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CHAPTER 1 
INTRODUCTION 



The General Theory of Relativity, developed by Albert Einstein at the beginning 
of the twentieth century, has been successfully used by scientists for many years to 
explain physical phenomena and to make predictions about physical systems. One of 
the most exciting predictions of General Relativity is the existence of gravitational ra- 
diation, which can be thought of as a wave-like distortion of spacetime. It is expected 
that gravitational radiation from remote astrophysical systems will yield very impor- 
tant information about those systems, thus answering many questions that scientists 
currently have. Even though the detection of gravitational waves has proved to be 
a challenging task, recent developments in technology have made scientists confident 
that gravitational waves will indeed be detected in the near future. Both earth-based 
detectors (such as LIGO and VIRGO) and space-based detectors (such as LISA) are 
expected to start operating in the next decade or two. For the data collected by 
these detectors to be useful, accurate information about the sources of gravitational 
radiation is required. Only if scientists know what the gravitational waves emitted by 
specific systems should look like, can they compare them with the patterns observed 
and draw conclusions about the physical characteristics of those systems. That is why 
there has been a lot of interest lately in predicting the gravitational radiation emitted 
by different astrophysical systems. 

One astrophysical system for which gravitational waves are expected to be de- 
tected, specifically by space-based detectors, is the binary system of a small neutron 
star or black hole (of mass equal to a few times the mass of the Sun) and a super- 
massive black hole (of mass equal to a few million times the mass of the Sun) [1]. 
Supermassive black holes are believed to exist at the centers of many galaxies, in- 
cluding our own. It is also believed that their strong gravitational field can capture 
smaller stars and black holes, which then move toward the supermassive black hole, 
until they are absorbed by it. The exact evolution of the system with time can give 
the pattern of the gravitational radiation that the system is expected to emit. Know- 
ing that pattern is crucial in distinguishing whether the gravitational waves measured 
by the detector come from such a system. It can also help extract information about 
many characteristics of the system, such as the masses and the angular momenta of 
the two components of the binary system. 

To determine how the system evolves, the path of the smaller star or black 
hole needs to be calculated. That path is largely determined by the influence of the 
gravitational field of the supermassive black hole on the small star or black hole; it is 
also affected by the interaction of the small star or black hole with its own gravitational 
field. This latter interaction, commonly referred to as the self-force, is responsible for 
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the radiation reaction effects, which cause the decay of the small star's or black hole's 
orbit toward the supermassive black hole. 

Even though my motivation for studying the self-force effects stems from the 
need to know the evolution of this system, radiation reaction effects are present in 
other systems that are easier to deal with mathematically. One can think, for example, 
of a particle of a certain scalar charge, which moves in a Schwarzschild background 
spacetime and creates its own scalar field. There is also the case of a particle that 
carries an electric charge and creates an electromagnetic field as it moves in spacetime. 
The scalar field of the particle in the first case and the electromagnetic field of the 
electric charge in the second case will affect the motion of each particle, causing its 
worldline to differ from what it would be if radiation reaction effects were not present. 
It is useful to predict the evolution of such systems, mainly because they can give an 
idea of how the more difficult system can be handled successfully and not because 
they are reahstic systems, since they are not expected to be observed in nature. 

In the past, various efforts have been made to develop a concrete scheme for 
including the radiation reaction effects for scalar, electromagnetic and gravitational 
fields in the equations of motion. One approach that has been used extensively in- 
volves calculating the fiux, both at infinity and at the event horizon of the central 
supermassive black hole, of quantities (such as the energy and momentum of the par- 
ticle) that are constants of the motion when radiation reaction is not present [2-9]. 
That flux is then associated with the rate of change of those quantities at the location 
of the particle and the evolution of the particle's motion can be predicted. However, 
this approach is generally not applicable in the case in which the central black hole 
is a rotating black hole, except under some special circumstances [10]. One reason is 
that the evolution of the Carter constant [11], which is a constant of the motion that 
is not derived by a Killing vector, cannot be calculated using the flux at inflnity and 
at the event horizon. In addition, this approach does not take into account the very 
significant non-dissipative effects of the self-force that the particle's field exerts on the 
particle [12]. That is why it is important for the self-force, and not just its radiation 
reaction effects, to be calculated. The calculation of the self-force is the subject of 
this dissertation. 

In Chapters 2 and 3 of this dissertation, I present the innovative work of Dirac 
[13] for the radiation reaction effects on an electron moving in flat spacetime and 
the theoretical generalization of it to the case of curved spacetime by DeWitt and 
Brehme [14]. In Chapter 4, I present different methods for practical calculations of 
the self-force that have been proposed in the past, with emphasis on one in particular. 
In Chapters 5-8 of this dissertation, I present an implementation of that particular 
method for calculating the radiation reaction effects for a scalar, an electromagnetic 
and a gravitational field. 



CHAPTER 2 
RADIATION REACTION IN FLAT SPACETIME 



The first successful attempt to provide an expression for tiie radiation reaction 
effects on a particle in special relativity was made by Dirac in 1938 [13]. In his famous 
paper, he studied the motion of an electron moving in flat spacetime by using the 
concept of conservation of energy and momentum along the electron's worldline. In 
this chapter 1 present the decomposition of the electromagnetic field and the derivation 
of the equations of motion, which include the radiation damping effects on the electron. 
The scheme for calculating the radiation reaction effects in curved spacetime presented 
in this dissertation extends Dirac's work. 

For the equations to look as simple as possible, the speed of light is set equal to 
1 in this chapter. 

It is noted that the metric signature used by Dirac is (H ) and that is what 

is used in this chapter, in order to keep the results identical to those derived by Dirac. 
In all subsequent chapters, however, the signature is changed to ( — h ++), in order 
to adhere to the currently widely used convention. Thus, some of the equations given 
in this chapter change in subsequent chapters. Since none of the results described in 
this chapter are used in the chapters that follow, that should not cause any confusion 
at all. 

2.1 Definitions of tlie Fields 

Dirac's analysis starts by assuming that an electron of electric charge q is moving 
in fiat spacetime, so that the background metric is the Minkowski metric 77^;^. The 
worldline of the electron is denoted by 

= ^^(«) (2-1) 

where s is the proper time along the worldline. The electromagnetic 4- vector potential 
that the electron creates is assumed to satisfy the Lorentz gauge condition 

VA^ = (2.2) 

and obeys Maxwell's equations 

V'A^ = 4ttj^. (2.3) 

It is well-known that the Lorentz gauge condition leaves some arbitrariness for the 
electromagnetic potential. In Equation (2.3), is the charge-current density vector 
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of the electron, namely 
f dz 

jix^q -j^ S{xq - zo) S{xi - zi) 5{x2 - Z2) S{x3 - zs) ds. (2.4) 

The electromagnetic field associated with a 4-vector potential is, in general, given 
by the equation 

F^" = V'A'' - V^A". (2.5) 

It is clear that, unless the boundary conditions of the problem are specified. 
Equations (2.2) and (2.3) do not have a unique solution. In fact, adding any solution 
of the homogeneous equation 

V^A^ - (2.6) 

(which represents a source-free radiation field and obeys the Lorcntz gauge condition) 
to a solution of Equation (2.3), gives a different solution of Equation (2.3) that also 
obeys Equation (2.2). 

One solution of interest for the problem at hand is the retarded electromagnetic 
field F^^ created by the electron, which is derived by the well-known Lienard-Wiechert 
potentials Afg^ that obey Equations (2.2) and (2.3). Assuming that there is also a 
radiation field F^^ incident on the electron, the actual field in the neighborhood of 
the electron is the sum of those two fields 

= FZ + Fr. (2.7) 

Another interesting solution is the one given by the advanced potentials A'^^. It is rea- 
sonable to expect that the advanced potentials play a role symmetrical to the retarded 
potentials in this problem. If the outgoing radiation field leaving the neighborhood 
of the electron is F^'^, then 

F^:t = FZ + FZ- (2.8) 

The radiation emitted by the electron is given by the difference of the outgoing and 
the incoming radiation 

Using Equations (2.7) and (2.8) for the actual field, an alternative expression for the 
radiation field can be obtained 

FZ = FZ - FZ- (2.10) 

The important implication of Equation (2.10) is that, because the retarded and the 
advanced fields are uniquely determined by the worldline of the electron, so is the 
radiation emitted by the electron. It becomes clear from the analysis of Section (2.2) 
that the radiation field is responsible for the deviation of the electron's worldline from 
a background geodesic, because that field shows up in the equations of motion, in the 
radiation damping term. 

It is significant to reahze that Equation (2.10) for the radiation field is consistent 
with what one would expect for the radiation of an accelerating electron. To under- 
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stand that, one should recall that the radiation produced by an accelerating electron 

can be calculated using the retarded field Fl^^'^ that the electron creates, at very large 
distances away from the electron and at much later times than the time when the 
acceleration takes place. At those distances and times the advanced field vanishes 
and Equation (2.10) gives the correct result. The advantage of Equation (2.10) is 
that it gives the radiation produced by the electron at any point in spacetime, so it 
can be used for the radiation in the neighborhood of the electron. 

A field that is needed for the analysis of Section (2.2) is the difference of the 
average of the retarded and advanced fields from the actual field, denoted by f^'' 

r-F^t-liFZ + F^.)- (2.11) 
Using Equations (2.7) and (2.8) with Equation (2.11) 

r'' = liF,7 + FZ)- (2.12) 

Because the incoming and outgoing radiation fields are derived from potentials that 
satisfy the homogeneous Equation (2.6), the field f^" is sourceless. It is also free from 
singularities at the worldline of the electron. 

2.2 Derivation of the Equations of Motion 

The main concept used by Dirac in deriving the equations of motion for the 
electron moving in flat spacetime is the conservation of energy and momentum. The 
first step is to surround the worldline of the electron by a thin cylindrical tube of 
constant radius. This radius is assumed to be very small, specifically smaller than 
any length of physical significance in the problem. The objective is to calculate the 
fiow of energy and momentum across the three-dimensional surface of the tube, using 
the stress-energy tensor T^p calculated from the actual electromagnetic field. The 
stress-energy tensor is given by 

4 TT T;, = F^t ,u Fact p + ^ Fa^ F^^ (2.13) 

The fiow of energy and momentum out of the surface of the tube is equal to the 
difference in the energy and momentum at the two ends of the tube. 

In the following, dots over quantities denote differentiation with respect to the 
proper time s. To simplify the notation, (dz/j^/ds) is set equal to V/j,. One then obtains 
the equations 

i^X = l' (2.14) 
v^v^^ = 0, (2.15) 
Vfj^v^" + Vf,vi^ ^ 0. (2.16) 

To calculate the stress-energy tensor, the electromagnetic 4- vector potential and 
the electromagnetic field need to be calculated first. The retarded potential at a point 
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X/j, generated by an electron moving on the worldline Zi^{s) is given by 



"-ret /i 



(2.17) 



Z„ [X'' - Z") 

calculated at the retarded proper time, which is the value of s that solves the equation 

{x,-z,){x'' -z")^^. (2.18) 

An equivalent expression for the retarded potential, obtained from Equation (2.17) 
by using elementary properties of the 5 function, is 



/Tint 
Zu, 5[(x^ - Zi,) {x" - z")] 
-oo 



ds 



(2.19) 



where Tint is a proper time between the retarded and advanced proper times. By 
differentiating Equation (2.19) and using Equation (2.5), the retarded electromagnetic 
field is obtained 



ret jxv 



r-* d 

J-oo 



Zv y^' ii Zp^) 



b\{xx - zx) {x^ - z^)\ ds. (2.20) 



By using again some of the elementary properties of the 5 function, an equivalent 
expression for the retarded field is derived 



ret f^u 



d 



z\ [x^ — z^) ds 



Zu Zy) Zn (Xp^ Zn 



Zk {x'^ — z'^) 



{2.21) 



Again, all quantities are calculated at the retarded proper time. 

Since the goal is to calculate the stress-energy tensor at the two ends of the 
world-tube surrounding the worldline of the electron, it can be assumed that the 
point Xp is very close to the worldline. Specifically 



Xp = Zp{So) +-fp 



(2.22) 



where the are very small. Then, the fields can be Taylor-expanded in the 7^'s. 
In the expansions that follow, all the coefficients are taken at the proper time Sq. It 
can be assumed that the retarded proper time is Sq — a, where it is reasonable that a 
is a small positive quantity of the same order as j^. 

The detailed calculations needed to obtain the Taylor expansions of the electro- 
magnetic fields are tedious but straightforward. They can be found in Dirac's paper; 
here 1 present only the results of that calculation. 

Since is a space-like vector, it can be assumed that 7^7^ = — where e is 
a positive number. Taylor-expanding the retarded field and keeping only the terms 
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that do not vanish in the hmit £ — > gives 
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(2.23) 



The advanced field can be obtained by changing e to —e and changing the sign of the 
whole expression. By using the retarded and advanced fields in Equation (2.10) the 
radiation field on the worldline is obtained 



4 q 



(2.24) 



The actual field then becomes 
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+ 2^ ^{'•j'lJ.lv-Vvli,) 



(2.25) 



To calculate the energy and momentum flow out of the tube, the stress-energy 
tensor needs to be calculated. In fact, only its component along the direction of 7^ is 
necessary. Substituting Equation (2.25) into Equation (2.13): 



4 TT T^p 7^ — i^act ixv Fact ^ -^act o:/3 -^act 7/i 



^2 ^1 



/I _4 1 _2' '\\ — 2/1 ^ '\\' 

(2 ^ + 2 ^ ^^^^ ^'^ ~ 2 ^ ^ 2 ^^^^ 



(2.26) 



The flow of energy and momentum out of the surface of the tube is given by 
integrating this component of the stress-energy tensor over the surface of the tube. 
The result is 

\ee-^v^ - qv^U,)ds (2.27) 



/ 



where the integration is over the length of the tube. 

As mentioned earlier, this flow of energy-momentum depends only on the con- 
ditions at the two ends of the tube. That means that the integrand must be a perfect 
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differential, namely 

^q^s-^v^ - qv^f^^^B^ (2.28) 

for some 5^. Equation (2.15) and the fact that f^^ is antisymmetric in its indices put 
a restriction on B^, specifically 

v^Bf^ ^^q^ vi^ - qv^v, ^ = 0. (2.29) 

Thus, the simplest acceptable expression for B^ is 

for some k independent of s. From Equation (2.28) one obtains 

k^]^q^ -m (2.31) 

where m must be a constant independent of e in order for Equations (2.28) and (2.29) 
to have a well-defined behavior in the limit £ — > 0. Finally, the equations of motion 
for the electron are 

mv^ = qvyf^ (2.32) 

and m plays the part of the rest mass of the electron. This result is used in Section 
(2.3) to derive an expression for the radiation damping effects on the moving electron. 

2.3 Radiation Reaction 

The field defined in Section (2.1) allows the equations of motion (2.32) to 

be expressed in a simple form. However, in practical applications one would prefer 
to have the incident radiation field Ff^^ in the equations of motion, as F-^^^ is usually 
given. By substituting Equations (2.7) and (2.10) into Equation (2.11), an expression 
for f^^ is obtained, that involves the incident field and the radiation field: 

r -Ft:-r\ FZ- (2.33) 

Then, the equations of motion become 

mv, = qv,F,^^-' + ^qv, F,,,/. (2.34) 

The first term of the right-hand side of Equation (2.34) involves the incident 
radiation field and gives the work done by that field on the electron. The second 
term of the right-hand side involves the radiation field emitted by the electron and 
is the term of particular interest when one considers radiation reaction. This term 
gives the effect of the electron radiation on itself and is present even if there is no 
external radiation field, namely even if F-^ — 0. The fact that the radiation reaction 
effects on the electron can be described by using the electromagnetic field F^ is 



(2.30) 
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one of the important results of Dirac's work and is the idea used when motion of 
particles in curved spacetime is considered, to calculate radiation reaction effects on 
those particles. 

It is worth writing the equations of motion in terms of the characteristics of the 
worldline of the electron and the external radiation field, which can be achieved by 
using Equation (2.24) for the radiation field of the electron 

mv^^qv^F^^/ +'^q'^v^ + ^ v^v" v^. (2.35) 

It is an important and a very interesting feature of Equation (2.35) that, in addition 
to the first derivative of t)^, its second derivative shows up as well. A discussion of 
that fact and of some of its implications is presented in [13], but that is beyond the 
scope of this dissertation. 



CHAPTER 3 

RADIATION REACTION IN CURVED SPACETIME 



Dirac's study of the radiation reaction effects on an electron gave a relatively 
simple result, because the electron was assumed to be moving in flat spacetime. The 
analysis becomes significantly more complicated for particles moving in curved space- 
time. DeWitt and Brehme [14] were the first to study the motion of an electrically 
charged particle in a curved spacetime using Green's functions. Later, Mino, Sasaki 
and Tanaka [15] generalized that analysis for the case of a particle creating a grav- 
itational field while moving in curved spacetime. Since Green's functions are used 
extensively in subsequent chapters to determine the self-force effects on different par- 
ticles, it is appropriate to present the analyses of DeWitt and Brehme and of Mino, 
Sasaki and Tanaka before proceeding. 

In this and all subsequent chapters, geometrized units are used, meaning that 
Newton's gravitational constant and the speed of light are both set equal to 1. Note 
also that the signature of the metric is changed to ( — h -|— 1-) from now on, so minor 
changes in some of the equations previously mentioned should not be surprising. 



One of the main characteristics of any Green's function is that it connects two 
points in spacetime. It does that by propagating the effect of the source, from the point 
where that source is located (source point) to the point where the field needs to be 
calculated (field point). Since Green's functions are inherently nonlocal quantities, the 
discussion about them can be facilitated if some more elementary nonlocal quantities 
are introduced first. 

The most general class of nonlocal quantities is the class of bitensors, which 
are simply tensors whose indices refer to two points in spacetime. In the following, 
X denotes the field point and z denotes the source point. In order for the indices for 
each point to be easily identifiable, all unprimed indices refer to the field point x and 
all primed indices refer to the source point z. 

A very important quantity for the study of nonlocal properties of spacetime is 
the biscalar of geodesic interval s{x,z). It is the magnitude of the invariant distance 
between the points x and z as measured along a geodesic that joins them and is a 
non-negative quantity. It is defined by the equations 



3.1 Nonlocal Quantities 



V"'sVa'S = ±l, 



(3.1) 



lim s = 0, 



(3.2) 
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where V" denotes covariant differentiation with respect to the background metric 
g'^''. It is obvious that s must be symmetric under interchange of its two arguments, 
namely 

s{x, z) = s{z, x). (3.3) 

With the signature of the metric being ( — h ++), the interval s is spacelike when the 
+ sign holds in Equation (3.1) and timelike when the — sign holds. The points x for 
which s — define the null cone of z. 

It is more convenient to use a different quantity to measure the invariant distance 
between the source point and the field point, a quantity that is, however, related to 
s{x,z). That quantity is Synge's [16] world function and is defined by 

a{x,z)= ±^s\x,z). (3.4) 

By using the defining equations for s one can deduce that a has the properties 

^ VV V«(7 = ^ V"V Va'a = a, (3.5) 

lim(7 = 0. (3.6) 

Also, a is positive for spacelike intervals and negative for timelike intervals. 

Another very significant nonlocal quantity is the bivector of geodesic parallel 
displacement, denoted by gaa'{x, z). The defining equations for it are 

V^WV;3(7 = 0, V^'wV/3/(7 = (3.7) 

lim^/ = C'- (3.8) 

Equation (3.7) signifies that the covariant derivatives of cjaa' ai'e equal to zero in the 
directions tangent to the geodesic joining x and z. Equation (3.8) expresses the fact 
that Qaa' is equal to the Kronecker-5 when x — z. The bivector of geodesic parallel 
displacement also has the property that 

9aa'{x,z) = ga'a{z,x). (3.9) 

In the following, the determinant of gaa' is denoted by 

9 = - \ gaa' I, (3.10) 

and 

^HTa'l- (3.11) 

The effect of applying the bivector of geodesic parallel displacement to a local vector 
A° at z is a parallel transport of that vector from point z to point x, along the 
geodesic that connects the two points. The result is a local vector A" at point x. In 
general, gaa' can be used to transform any bitensor that has indices refering to the 
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two points X and 2; to a tensor whose indices refer to one point only, for example 

-l3 rpa' P' _ rpa (3 /o ir,\ 

9 a' 9 f}' 7 5e — 7 <5e- l^^-iZj 

How useful the bivector of geodesic parallel displacement is in studying nonlocal pror- 
erties of spacetime can be understood if expansions of a bitensor about one point are 
considered. In order for a bitensor to be expanded about a certain point in spacetime, 
all its indices must refer to that specific point. For bitensors for which that is not the 
case, their indices must first be homogenized by applying g^a' and then the expansion 
about the specific point can be taken. 

A bivector that is very useful for the Hadamard [17] expansion of the Green's 
functions is defined by 

Daa'ix,z) = —Va'\^a'0'ix,z) (3.13) 

and a biscalar relating to it is its determinant 

D^-\ D^^, I . (3.14) 

DeWitt and Brehme proved that 

lim D^^,{x, z) = gaa'{z) (3.15) 

which shows that the biscalar D is nonvanishing, at least when x and z are close to 
each other. In fact, D is the Jacobian of the transformation from the set of variables 
{z" , x"^}, which specify the geodesic between x and z in terms of its two end points, 
to the set of variables {2;" , V" cr}, which specify the geodesic in terms of one of its 
end points and the tangent to the geodesic at that end point. DeWitt and Brehme 
also showed that D obeys the differential equation 

L>-^ Va(L> W) =4. (3.16) 

Instead of the biscalar D, a different biscalar is used in the Hadamard expansions 
of the Green's functions. That biscalar is denoted by A and is defined as 

A{x,z)^[g{x,z)]-^D{x,z). (3.17) 

DeWitt and Brehme proved that A can be expanded in terms of derivatives of a and 
that the expansion is 

A = 1 + ^ i?"^ VaCr Vpa + O(s^), loYx^z (3.18) 

where R"^ is the background Ricci tensor. 

3.2 Green's Functions 

The goal of the work of DeWitt and Brehme was to study the radiation damping 
effects on a particle of a given electric charge that moves in curved spacetime. To 
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do that it is necessary to study the vector field that represents the electromagnetic 
potential created by the particle. However, the equations for the scalar field are 
less complicated and thus easier to deal with. Also, the lack of indices makes the 
results more transparent. For those reasons the scalar Green's functions and their 
corresponding fields are discussed first. 

3.2.1 Green's Functions for a ScalEir Field 

A point particle of scalar charge q which is moving on a worldline F : ^"'(r), 
where r is the proper time along the geodesic, creates a scalar field -0 which obeys 
Poisson's equation: 

V^V' = -4 7r^. (3.19) 
In Equation (3.19), q is the source function for the point particle. Specifically 

Q{y) = q J i-g)-'' - z{t)) dr (3.20) 

where y is some point in spacetime. It is desired to express the different solutions of 
Equation (3.19) as integrals containing Green's functions, namely 



ip{x)=4:T: / {-g)^ G{x,y) g{y)dy^ 

J spacetime (3 21) 

— Airq Jg[x,z{t)] dr, 



so the properties of the various Green's functions are emphasized in this section. 

One function of importance is the symmetric Green's function, G®^™(x, z), which 
satisfies the inhomogeneous differential equation 

W'^G'^'^ix, z) = -{-gy^^ S\x - z). (3.22) 

The Hadamard form [17] of this function is 

G'^'^ix, z) = — [U{x, z) d{a) - V{x, z) e(-a)] (3.23) 
Stt 

where © is the step function that equals 1 if the argument is greater than zero and 
equals otherwise. U {x, z) and V{x, z) are biscalars that are free of singularities and 
symmetric under the interchange of x and z, namely 

U{x,z) = U{z,x) (3.24) 



and 

V{x, z) = V{z, x). 



(3.25) 
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They can be determined by expanding the solution of Equation (3.22) in powers of cr, 
in the vicinity of the geodesic F. The biscalar U {x, z) satisfies the differential equation 

U-^{x, z) V"U{x, z) = ^A"^(x, z) V"A(a;, z) (3.26) 
and the boundary condition 

lim C/(x, z) = 1 (3.27) 

and is given by 

C/(x,^) = [A(x,z)]5 

1 o (3-28) 

= 1 + — V^(7 + 0(5^), X ^ r 

where s is the proper distance from the point a; to F measured along the spatial 
geodesic which is orthogonal to F. The biscalar V{x^ z) satisfies the homogeneous 
differential equation 

V^V{x,z) = (3.29) 

and is given by 

V{x, z) = --^ R{z) + 0{s), x^r. (3.30) 

It is noteworthy that the symmetric Green's function vanishes for o" > 0, that is for 
spacelike separation of the points x and z. Also, because both U{x, z) and V{x, z) 
are symmetric under the interchange of x and 2;, so is G^^™: 

G'^y'"(x,z) = G^y°^(z,x). (3.31) 

If (7 = da{x, z{T))/dT, the scalar field ■?/'^^'™(x) associated with the symmetric 
Green's function is 



+ 

Tret 



' -00 -J Tadv 



and consists of two parts. The first part is the one that contains the biscalar U {x, z) 
and the 5-function 5{a) and is referred to as the direct part. This is the term that 
corresponds to cr = 0, which is the part of the field that comes from the retarded 
and advanced proper times (rj-ct and Tadv respectively), namely the proper times that 
correspond to the intersection of the geodesic F with the past and future null cone of 
the point x. In other words, the direct part has the same singularity on the null cone 
that the symmetric scalar field has in fiat spacetime. The second part is the one that 
contains the biscalar V{x, z) and the step-function ©(— cr) and is referred to as the 
tail part. This term gives the part of the field that corresponds to cr < 0, which is the 



15 



part that comes from the interior of the past and future null cone of x. It is the part 

of the field that is due to the curvature of spacctime and vanishes in flat spacetime. 

The symmetric Green's function can be separated into the retarded and ad- 
vanced parts, which constitute two very important Green's functions themselves. 
Specifically it can be written as 

G^^'^ix, z)^\ [G'^\x, z) + G^^{x, z)] (3.33) 

where the retarded and advanced Green's functions are given by 

G"^"*(x, ^) = 2 e[E(x), z\ G^'^'^ix, z), (3.34) 
G^dv^^^ _ 2 Q[z, E{x)] G'^^^ix, z). (3.35) 

In Equations (3.34) and (3.35), S(a;) is any spacelike hypersurface that contains x. 
The step-function 0[E(a;),2;] = 1 — 0[2;, E(a;)] is equal to 1 when z lies to the past 
of S(a;) and vanishes when z fies to the future of S(a;). Both the retarded and the 
advanced Green's functions satisfy the inhomogeneous differential equation 

V'G'""*(a;,^) = \j''G'^\x,z) = -{-g)''^ 5\x,z). (3.36) 

They also have the property that 

G''\x,z) = G'"^''{z,x). (3.37) 

The retarded and advanced scalar fields, both solutions of the inhomogeneous 
differential Equation (3.19), can also be expressed as integrals of the respective Green's 
functions 

^P''\x) ^ Anq [ G''''[x,z{T)]dT = 2[i^^^^l -q nW{x,z)dT, (3.38) 

Jr I 2a i Tret J-oo 

^^-{x) = 4 7rq f G'^''[x,z{T)]dT = 2[^-^^^^l -q r"v{x,z)dT. (3.39) 

Jr L 2(7 JTadv iTadv 

The retarded field V'''''* (x) is the actual scalar field that results from the scalar particle 
and is singular at the location of the particle. 

The last scalar Green's function that is significant for the analysis that follows 
is the radiative Green's function G^^'^, which is defined in terms of the retarded and 
advanced Green's functions in a manner analogous to the radiation field defined by 
Dirac in fiat spacetime, namely 



G'^^'^ix, z) = G'^\x, z) - G'^^ix, z). (3.40) 

It is obvious by this definition that G^^{x, z) satisfies the homogeneous differential 
equation 

W'^G""^{x,z)^0. (3.41) 
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By Equation (3.37) one can infer that G^^'^{x, z) is antisymmetric under the inter- 
change of its arguments 

G''"^{x,z) = -G''''^{z,x). (3.42) 

The corresponding radiation field equals the difference between the retarded and 
advanced fields 

^''^{x) = ^l^''\x) - ^"^""{x) ^4:nq J G'^[x, z{t)] dr (3.43) 

and is a sourceless field because it satisfies the homogeneous differential equation 

V^^"^{x) = 0. (3.44) 

Consequently, ip'^'^'^ is smooth and differentiable everywhere in space. 

The Green's functions and the corresponding scalar fields that were mentioned 
in this section are used extensively in Chapter 4, where the self-force on scalar particles 
is discussed. 

3.2.2 Green's Functions for a Vector Field 

Now that the scalar Green's functions and their properties have been discussed, 

it is straightforward to define analogous Green's functions for the electromagnetic 
vector potential A" generated by an electrically charged particle moving on a geodesic 
r : z"' (r). Assuming that A" satisfies the Lorentz gauge condition 

Va^" = 0, (3.45) 

Maxwell's equations for it become 

V^^° - R'^pA^ = -4 TT J°, (3.46) 

and, again, the goal is to express the vector potential as an integral that contains a 
Green's function and the source J", specifically 

A,(x)=47r [ {-gY^G^p{x,y)jP{y)d^y. (3.47) 

J spacetime 

The symmetric Green's function for a vector field obeys the equation 

V^Cr z) - R^Gj^ix, z)^-g--^ w z) (3.48) 
and the Hadamard form of it is 

GTc^i^^z) = ^[Uaa'{x,z) 5{a) - V,,^,{x,z) e(-a)]. (3.49) 
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The biscalar Uaa' is given by the differential equation 

(2 V^Uaa' + Uaa' A"^ V^A) V pa = (3.50) 
with the boundary condition 

limUaa'{x,z) = gaa'{z). (3.51) 

x—*z 

As shown in [14], the solution is unique 

Uaa'{x,z) = ^/A{X,Z) Qaa'ix^z) 



[1 + ^ R"'^' Vya + 0{s')] gaa', for X 



(3.52) 



where Equation (3.18) was used to derive the final expression. For the bi vector Vaa' 
DeWitt and Brehme prove that 

\im Vaa'ix, ^) = ^ 9a ' i^a' p' " ^ 9a' p' R) + 0(s), for x^V. (3.53) 
The symmetric electromagnetic vector potential is calculated by 

= IJ {-9f^{[A{x,yt^ -gap{x,y) 6{a) - Vap{x,y) e(-a)} J^{y) d'y. 

(3.54) 

Just as in the case of the symmetric scalar field, the symmetric electromagnetic poten- 
tial also consists of the direct part, that contains the 5-function and comes from the 
retarded and advanced proper times, and the tail part, that contains the 0-function 
and is the contribution from within the past and future null cone. An interesting fea- 
ture is the appearance of the bivector of geodesic parallel displacement in the direct 
part. It signifies that the electromagnetic radiation is parallel propagated along the 
null geodesic that connects the points x and z. 

The symmetric Green's function can again be separated into the retarded and 
advanced parts 

^) = \ [GfA^. z) + Gt:\x, z)] (3.55) 
which are given by the equations 

Gf^x, z) = 2 e[E(x), z\ G^^{x, z) (3.56) 

Gt:,{x, z) = 2 Q[z, E(a:)] G^'ix, z). (3.57) 
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They are solutions of the inhomogeneous differential equation 

W^Gf^,{x,z)-R^^GfAx,z) = V'G^:,{x,z)-Rj^Gf:,{x,z) = --g-^^^,{x,z)5\x,z) 

(3.58) 

and have the property that 

GTA^,z) = G^l{z,x). (3.59) 
The retarded and advanced electromagnetic potentials are then given by 

Af{x)^An I {-gY^G%{x,y)J^{y)d'y (3.60) 

J spacetime 

A^^ix)^An I {-g)'^G^^{x,y)jP{y)d'y (3.61) 

J spacetime 

and A^^^{x) is the actual potential generated by the moving charged particle. 

The radiative Green's function is defined in terms of the retarded and advanced 
Green's functions by the equation 

G-^,(x, z) = G-*,(x, z) - Gf:,{x, z) (3.62) 

and satisfies the homogeneous differential equation 

V^G-^,(x,^) - K^GX\x,z) =0. (3.63) 

Using the symmetry expressed in Equation (3.59) for the retarded and advanced 
Green's functions it can be inferred that 

G^:^,{x,z) = -G^:;i{z,x). (3.64) 

The corresponding radiation vector potential is equal to the difference between 
the retarded and advanced vector potentials 

A^:'{x) = Afix) - A^\x) 

= 4 7r / {-gy^G^:^{x,y)J^{y)d'y ^^"^^^ 

J spacetime 

and is a source-free field, because it satisfies the homogeneous differential equation 

V''A'^'^{x) - R£A'p^{x) = 0. (3.66) 

The important quantities for the motion of a charged particle are not the vec- 
tor potentials but rather the electromagnetic fields. The electromagnetic field that 
corresponds to a vector potential A^ is calculated by: 

Fafsix) = VpA^ix) - VaAf^ix). (3.67) 
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For completeness, the definitions of tlie symmetric and radiation fields are given, in 
terms of the retarded and advanced fields: 



er = ^(^a| + (3-68) 



Kf = ^a? - ^a?- (3.69) 



3.2.3 Green's Functions for a Gravitational Field 

Even though the work of DeWitt and Brehme did not cover the motion of 
a massive particle creating its own gravitational field, the Green's functions for a 

gravitational field arc presented here, for completeness. A more detailed description 
is given by Mino, Sasaki and Tanaka in [15] and it is that analysis that is followed 
here. 

Assuming that a massive pointlike particle is moving on a geodesic F : (r) 
and is inducing a perturbation h^uix) on the background, the trace- reversed metric 
perturbation is 

hu,^(x) = /i^^(x) - -g^,y{x)h{x) (3.70) 

where h{x) and h{x) are the traces of h^^{x) and of h^^{x) respectively. It is assumed 
that the trace-reversed metric perturbation obeys the harmonic gauge condition 

Wji^^ix) = 0. (3.71) 

In this gauge, the linearized Einstein equations [18] become: 

V + 2i?j/ h^s = -IGttT,^ (3.72) 

to first order in the metric perturbation. The Green's functions of interest in this 
case are bitensors and are used to express the trace-reversed metric perturbation in 
integral form 

h'^^ix) = 167r I {-g)'^G'^^^\x,y) T,s{y) d'y. (3.73) 

J spacetime 

The symmetric Green's function G"^^'^' {x, z) satisfies the differential equation 
V^G^/^'^'(x, z) + 2 i?"/,(x) G^;2:''i^: = -2 r'^"{x, z) f^'{x, z) (3.74) 

and the Hadamard form for it is 

G"fj'^'(x, z) = ^[C/"^^'^'(x, z)5{a) - V^^^'^'ix, z)Q{-a)]. (3.75) 
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The bitensor U°'^'^'^'{x, z) is the solution of the homogeneous differential equation 

[2 V«C/"^^'^' (x, z) + ^ U""^^'^' {x, z)] V^a = (3.76) 
with the boundary condition 

lim W'^^'^'ix, z) = lim 2^^'("(x, z) f^^'{x, z). (3.77) 

x^z x^z 

Mino, Sasaki and Tanaka [15] prove that the solution is 

C/"^^'^'(x, z) = 2 ^^'("(x, z) f^^'ix, z) \/^{x,z). (3.78) 
The bitensor V^^^'^' (x, ^) is divergence-free 

V^y'*^^''^'(x,z) = (3.79) 
and satisfies the homogeneous differential equation 

y2ya/37'5' (^^ z) + 2 R""/^ V'"^'^' (x, z) = (3.80) 

with the boundary condition 

limy"^^''^'(x,z) = 0. (3.81) 

fT— >0 

It is proven in [15] that the solution is 

K/37'5'(a^, z) = -9a' ffi R^'a'5'0'{z) + 0{r), for X ^ r. (3.82) 
The symmetric gravitational field is obtained by 

h%^^{x) = IGtt I (-5)^ G,f/(a;, y) T,s{y) A 

J spacetime 

= 2 / {-gf^ [U'^f'^\x,y)6{a) - V'^^^\x,y)ei-a)] T,,{y) d'y. 

J spacetime 

Its direct part is the part that contains the 5-function and gives the contribution from 
the null cone of point x. Its tail part is the part that contains the 0-function and 
gives the contribution that comes from within the null cone of point x. 

The symmetric Green's function can be separated to the retarded and advanced 
Green's functions 

Gtf2'''(x, z) = \ [G-^C'^ix, z) + Gt:''\x, z)] (3.84) 
which are given by the equations 

Gt'''\x, z) = 2 e[E(x), z] Gl^X'\x, z) (3.85) 
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z) = 2 e[z, E(x)] G',fj'^'(x, z). (3.86) 
They are solutions of the inhomogeneous differential equation 

V^G^;f'^'(x, z) + 2 R^^ix) GZ'''\x, z) = 

(3.87) 

and have the property that 

G:i'^''\x,z) = G]Z''^iz,x). (3.88) 
The retarded and advanced gravitational fields are given by 

Kft = 167r / i-gf^ G^i-^'ix, y) T,s{y) d'y (3.89) 

J spacetime 

^1 = / (-9)'^ G'^ITi^^ y) TAV) d'y (3.90) 

J spacetime 

and is the actual trace-reversed gravitational perturbation induced by the moving 
particle. 

The radiative Green's function is defined in terms of the retarded and advanced 
Green's functions by the equation 

gZ''' ^) = Grt" ^) - g:!:'' {x, z) (3.91) 

and satisfies the homogeneous differential equation 

V'G^^''\x, z) + 2 R\^M G^J'ix, z) = 0. (3.92) 

Using the symmetry expressed in Equation (3.88) for the retarded and advanced 
Green's functions, a symmetry property for the radiative Green's function can be 
derived, namely 

G:^'''\x,z) = -Gt^-^iz,x). (3.93) 

Finally, the radiative gravitational perturbation is given by the difference be- 
tween the retarded and advanced gravitational perturbations 



= 167r / (-^)iG'f/(x,y)r,,(y)A 

J spacetime 



(3.94) 



' spacetime 

and is a solution of the homogeneous differential equation 



V'K^,{x) + 2 i?"/, hZix) = 0. (3.95) 
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3.3 Motion of a Charged Particle in Curved Spacetime 

Following DeWitt and Brehme's analysis, I present the equations of motion and 
the radiation reaction effects for a particle carrying an electric charge e. Studying the 
vector case is useful in comparing the result derived in curved spacetime with Dirac's 
result for flat spacetime. Similar analyses can be done for the scalar case and the 
gravitational case. Those analyses can be found in a recent review of the subject by 
Poisson [19]. 

3.3.1 Equations of Motion 

Let's assume that a particle of bare mass mo and electric charge e is moving in 
a curved spacetime of metric g^^ on a worldline described by T : z{t), where r is the 
proper time. 

In the following, the notation u"'' — {dz°'' /dr) = z""' is used. Also, the symbol 
(D/dr) is used to denote absolute covariant differentiation with respect to the proper 
time r along the worldline. 

The particle generates an electromagnetic vector potential A°'{x) and an elec- 
tromagnetic field F'^^lx) given by Equation (3.67). The action of the system contains 
three terms; the first term comes from the particle, the second from the electromag- 
netic field and the third from the interaction between the two: 

S Sparticle ~l~ S field ~l~ Sinteractiom (3.96) 

where 

Sparticle ^ -mo J [-ga'i3'{z) u"^' u'^']^ dr (3.97) 

S field = - J^ap d^'x (3.98) 

S..ter.ct.o.-e jAAz)u-'dr. (3.99) 

In determining the equations of motion, the stationary action principle can be 
applied to S for variations in the potential A"' and for variations in the worldline z°^ 
separately. Varying A"" gives 

VflF"^ = 47rJ" (3.100) 
where J" is the current density defined as 

J'^iv) = ejjy{y,z) u^' 5\y,z) dr. (3.101) 

The electromagnetic field Fq,^ is invariant under a gauge transformation of the form 
Aa ^ Aa-\- VaA(x), wherc A(x) is any scalar function of x. For convenience, A(x) 
can be chosen in such a way that the electromagnetic potential satisfies the Lorentz 
gauge condition 

Va^" = (3.102) 
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in which case Equation (3.100) for the electromagnetic field becomes 



V'^A'^ix) - R'lA^ix) = -47rJ"(a;). 



(3.103) 



Equation (3.103) allows one to determine the electromagnetic field, once the worldline 
of the charged particle is known. 

On the other hand, varying the worldline gives 



Equation (3.104) allows one to calculate the worldline of the particle, once the elec- 
tromagnetic field is specified. 

As is well-known, the electromagnetic field given by Equation (3.103) is diver- 
gent at the location of the particle and cannot be used in Equation (3.104). An 
alternative method for studying the motion of the particle is to take advantage of the 
conservation of energy and momentum for the particle and the field, the same method 
that was used by Dirac in his analysis of the motion of an electron in fiat spacetime. 

3.3.2 Calculation of the Fields 

It becomes obvious in Section (3.3.3) that in order for the conservation of energy 
and momentum to be implemented along the worldline of the charged particle, it 
is necessary for the retarded and advanced fields generated by the particle to be 
calculated. A brief description of the calculation is given here; the details can be 
found in [14]. 

In the following, Tj-^t and Tadv are the proper times along the worldline of the 
particle, where the past and future null cones of x respectively intersect that worldline. 
They are referred to as retarded and advanced proper times. Also, ts is the proper 
time at the point where the hypersurface S(x) intersects the worldline of the particle. 
All these points are shown in Figure (3.1). 



eF'^Uz) vP' . 



(3.104) 




Figure 3.1: Retarded and advanced proper times 
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Starting with Equation (3.60) for the retarded potential, and using the Hada- 
mard expansion for the Green's function, a simple expression for A^^^ can be obtained 

ATix) = 47: J d'y {-g)'^ G%{x,y) J^{y) 

^47re / A(-^)^2e[E(x),y] P<^^M K^) - V'^^i^^v) ^i'^)] ^ 

J OTT 

X J^dTgUy,z) 5\y,z) 
Vaa'{x, z) 6(— cr) ii" dr 

-oo 

Vaa' {x, z) m"' dr. 



^ ^lU UP J Tret 



— e 



(3.105) 



By a similar calculation, starting with Equation (3.61), the advanced potential can 
be calculated 



/•+00 

y 



Vaa'(2^, ^) dr. 



(3.106) 



Equations (3.105) and (3.106) give the covariant Lienard-Wiechert potentials. 

The retarded and advanced electromagnetic fields can be calculated by Equation 
(3.67) if the respective potentials are differentiated. The result is 



X 



dr 

X (Vy(7 M^')-^ 



Tret / Tadv 



±oo 



(V^Ka' - V.y^aO li'*' dr. 



Tret / Tadv 



(3.107) 



In Equation (3.107), the lower of the two signs corresponds to the retarded field and 
the upper sign corresponds to the advanced field. This is an important expression 
because it gives the electromagnetic field at point x in terms of the characteristics of 
the worldline of the charged particle. The significance of this becomes clear when the 
final equations of motion of the charged particle are written down. 

As becomes obvious in Section (3.3.3), when calculating the balance of energy 
and momentum between the particle and the field along the worldline, the values of 
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the fields very close to the worldline need to be known. Expanding the retarded and 
advanced fields in powers of a and its derivatives is a tedious but straightforward 
calculation, which yields a complicated expression for the two fields, Equation (5.12) 
of [14]. 

Following Dirac's ideas, the total field in the vicinity of the charged particle can 
be written as the sum of the retarded field created by the particle plus the external 
incoming field 

F^p = + Ftp (3.108) 

which, by using Equations (3.68) and (3.69), can be written in terms of the radiation 
and symmetric fields 

^a/J = i^"^ + ^F^/ + i^X- (3-109) 

This breaking up of the total field is very important because, after implementing the 
conservation of energy and momentum along the worldline of the particle, only the 
radiation field F^^ and the external incoming field F^^ show up in the equations of 
motion of the charged particle. 

Since the radiation field is the difference between the retarded and advanced 
fields, only two terms of Equation (3.107) are of interest here. Specifically, it is 
sufficient to write 



F',f^^ = e {g^rg^p. - g. 



'9 1^13' 



2 fD^a\~^ Du^' 



1 

± - 
2 



/itoo 



(3.110) 



+ F^"^ + (higher order terms) 



where, again, the lower sign corresponds to the retarded field and the upper sign 
corresponds to the advanced field. In Equation (3.110), the notation 

fafia'ix, Z) = VpVaa'ix, z) - VaVga/(a;, z) (3.111) 

is used. The term F^^'^^ is, in fact, a sum of terms and is the same for the retarded 
and the advanced fields, so it vanishes when the difference of the two is taken. 

3.3.3 Conservation of Energy and Momentum 

The first step, when considering the conservation of energy and momentum, is 
the construction of a cylindrical surface that surrounds the worldline of the charged 
particle. That surface is referred to as the world-tube from now on. As in Dirac's 
analysis, the radius e of the world-tube is very small, because the balance of energy 
and momentum between the particle and the field needs to be calculated in the vicinity 
of the particle's trajectory. DeWitt and Brehme [14] explain the procedure for the 
construction of the world-tube in great detail. For the needs of the present calculation, 
the world-tube can be thought of as the three-dimensional surface that is generated 
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by the particle if the particle is surrounded by a sphere of a very small radius and is 
then left to move along its worldline, starting at proper time Ti and ending at proper 
time T2. In the following, the cylindrical surface is denoted by E, its two end caps are 
denoted by Ei (corresponding to the earliest proper time ri) and E2 (corresponding 
to the latest proper time T2). The total volume enclosed by the world-tube is denoted 
by 14- 

The stress-energy tensor of the system consists of one part that comes from 
the particle and one part that comes from the field 

T'^P^Tf + Tf (3.112) 

where 

Tf ^moj6-2 g''^,g^^, u'^'u^^' 6^ dr (3.113) 

Tf = (F"F^^ - \ g'^'^F^sF^^) (3.114) 

4:71 ' 4 

where F^js is the total field in the vicinity of the particle. It is the sum of the field 
generated by the particle plus any external incoming field that is present. 

The conservation of energy and momentum for the system can be expressed as 

V^T"^ = 0. (3.115) 

When integrating Equation (3.115) over the world-tube, special attention should be 
given to the fact that the integral (J V pT^-^d^x) is not a vector in the curved spacetime 
that is being considered here. That means that Gauss's theorem cannot be used as 
usual, to convert the integral over the volume Vt to an integral over the surface of the 
world-tube. The natural way to overcome this difficulty is to consider the integral 
/ Qa''^ fl'°'^d'^x instead. In this integral, the bivector of geodesic parallel displacement 
is used, so that any contributions to the integral by the point x are referred back to 
the fixed point z, which is assumed to correspond to proper time r. That integral is, 
then, a contravariant vector at point so Gauss's theorem can be used. That gives 

( X + ^ + ^ ) d^P - J^^ (V/3^a"') T"-^ d'x = (3.116) 

where dE^ is the surface element along the surface of the world-tube. 

The limit of the radius of the world-tube going to zero is taken next. Also the 
integral over the surface E of the world-tube is expressed as a double integral, over 
the proper time and over the solid angle. In addition, the proper times ri and T2 are 
let to approach r. If their infinitesimal difference is denoted by dr. Equation (3.116) 
becomes 

mo——dT + lim / gfT'^'^dEa = 0. (3.117) 

It is in the calculation of the integral of the second term that the symmetric field and 
the radiation field mentioned in Equations (3.68), (3.69) and (3.109) are very useful. 
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As was mentioned earlier, the result contains only the incoming external field, 
the radiation field and the tail field and is 



m- 



eF'''''p,u^' + \er'^'''p,u^' + \e\^' f f%,y{z{T), z{t'))u^' {T')dT' (3.118) 
2 2 



dr 

where m is the renormalized observed mass of the charged particle 

m — rrio + \im-e^e~^ . (3.119) 

e— >0 2 

An equivalent expression is derived by using Equation (3.110) and the definition 
of the radiation field given in Equation (3.69), and it is 



+ \e\^' rp'A<r), z{t'))u^' {T')dT' . 



Equation (3.120) contains only the physical characteristics of the worldline of the 
charged particle and the incoming external field. 



3.3.4 Radiation Reaction 

A short discussion about Equations (3.118) and (3.120) is in order. First it is 
significant to notice that, except for the last term that contains the quantity fapa'i 
these equations are identical to Equations (2.34) and (2.35) derived by Dirac for fiat 
spacetimc. Since that term contains the biscalar V^a', it comes from the tail part of 
the retarded electromagnetic field. 

Also, it is interesting that the integral that shows up in these two equations 
needs to be calculated from — oo to the proper time r, meaning that knowledge of the 
entire past history of the particle is required. That is something that renders the use 
of Equation (3.120) impractical and is a point that is also discussed in the following 
chapters. 



CHAPTER 4 
SELF-FORCE 



In general, the self-force can be due to a scalar, an electromagnetic or a gravita- 
tional field. Two different metfiods for the self-force calculation are presented in this 
chapter. The analysis is given in detail for the scalar self-force, due to the simplicity 
of the notation. The results are also given for the electromagnetic and gravitational 
self-force. 

The first method, in which the direct and tail parts of the retarded fields are 
used, was described by DeWitt and Brehme in [14] for the scalar and electromagnetic 
fields and by Mino, Sasaki and Tanaka in [15] for the gravitational field. An axiomatic 
approach to this method was also presented by Quinn [20] for the scalar field and by 
Quinn and Wald [21] for the electromagnetic and gravitational fields. 

The second method, in which the singular-source part and the regular remainder 
of the retarded fields are used, was proposed by Detweiler and Whiting in [22], where 
they described it for the scalar, electromagnetic and gravitational fields. 

As in Chapter 3, z denotes a point along the worldline of the moving particle 
and X any point in spacetime. The limit of x ^ F is the coincidence limit. The 
retarded and advanced proper times. Tret and Tadv, are the proper times at the points 
where the null cone of x intersects the worldline F (see Figure (3.1)). Primed indices 
refer to z and unprimed ones refer to x. 



4.1 Scaletr Self- Force 

A particle of scalar charge q is assumed to be moving in a background spacetime, 
described by the background metric gah- For simplicity it can be assumed that there is 
no external scalar field and consequently, if the scalar charge q is small, the lowest order 
approximation for the worldline F : z"' (r) of the particle is a background geodesic, 
where r is the proper time. But the scalar field ^ generated by the particle interacts 
with the particle, inducing a self-force on it. The self-force gives F an acceleration, 
so, to higher order, the particle does not trace a background geodesic. 



4.1.1 Direct and Tail Fields 

The discussion of the scalar field self-force in [14] and [20] shows that the self- 
force on the particle can be calculated by 



= q 



a„;,self 



particle 



(4.1) 
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where ■0^^'^ is the scalar field that interacts with the particle and is given by 



Tret 



^-^(x) = - [^^1 - g / V{x, z{r)) dr. (4.2) 

L Za J Tret J_oo 

It can be inferred from the discussion of the scalar Green's functions of Section 
(3.2.1) that the first term of the right-hand side of Equation (4.2) comes from the 
direct part of the retarded field generated by the moving particle. This first term 
is finite and differentiable at the location of the particle and its contribution to the 
self-force is the curved-spacetime generalization of the Abraham-Lorentz-Dirac force 
of flat spacetime. This contribution results from the moving particle's acceleration, if 
the worldline F is not a geodesic, and can be expressed in terms of the acceleration of 
r and components of the Riemann tensor [14, 15, 20, 21]. 

The integral term of the right-hand side of Equation (4.2) comes from the tail 
part of the retarded field of the moving particle. Its contribution to the self-force 
represents the result of the scattering of the retarded field of the particle, due to the 
curvature of spacetime. Taking the derivative of the tail part of ip^^^^ gives one term 
that comes from the implicit dependence of the retarded proper time on x and one 
term that comes from the dependence of the integrand on x, as was calculated by 
Quinn in [20]: 

Va{-g r\{x,z{T))dT] 
J —oo 

/Trot 
{VaV{x,z{T))} dr 
-oo 



1 










Tret J — 



(4.3) 



qR{x}_ _ 
12 r ^""^ 



/Tret 
{VaV(x,z(r))}dr. 
-oo 



In Equation (4.3), r is the proper distance from x to F measured along the spatial 
geodesic that is orthogonal to F. 

The spatial part of the first term of the derivative calculated in Equation (4.3) 
is not well defined when x is on F, unless R{x) — there. That makes this approach 
problematic since, for the self-force to be calculated, the derivative must be calculated 
at the location of the particle. To overcome this difficulty, Va^^*^'^ can be averaged 
over a small spatial 2-sphere surrounding the particle, thus removing the spatial part. 
Then the limit of the radius of the 2-sphere going to zero can be taken and a finite 
contribution to the self-force is obtained. 

Another complication comes from the fact that, in order for the contribution of 
the tail term to the self-force to be calculated, knowledge of the entire past history of 
the moving particle is necessary, as is clear from the integral term in Equation (4.3). 
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It is also important to note that the field ^0 , just hke the field ijj^^ , is not a 
physically well-understood field. That is because it is not a solution to any particular 
differential equation. 

4.1.2 The S-Field and the R-Field 

In order to provide an alternative scalar field that gives the self-force on the 
moving particle, Detweiler and Whiting [22] took advantage of the fact that adding a 
homogeneous solution of a differential equation to an inhomogeneous solution of the 
same differential equation, gives a new inhomogeneous solution. Specifically, since 
the biscalar V{x,z) is a symmetric homogeneous solution of the scalar differential 
Equation (3.29), it can be added to the symmetric Green's function G^^^{x, z) to give 
a new symmetric inhomogeneous solution G^{x, z) of Equation (3.22) 

G^{x, z) = G^y"'{x, z) + ^V{x, z) 

OTT 

- i-[[/(x, z)d{a) - V{x, z)Qi-a) + Vix, z)] (4.4) 
= ^[U{x,z)S{a)+V{x,z)e{a)], 

OTT 

which obeys 

V^G^ix, z) = -{^g)U\x - z). (4.5) 

The superscript S is used to indicate that this Green's function obeys the differential 
equation that contains the Source term. This new symmetric Green's function has 
support on the null cone, coming from the term that contains the 5-function, and 
outside of the null cone, at spacelike separated points, coming from the term that 
contains ©(o"). It has no support within the null cone, as shown in Figure (4.1). 
The local expansion of the biscalar V{x,z) given in Equation (3.30) is sufficient for 
the purposes of this dissertation, since the singular Green's function is only used for 
points X close to the worldline of the particle. 

The Regular-Remainder Green's function is defined in terms of the retarded and 
singular Green's functions as 

G^{x, z) = G'^\x, z) - G^{x, z) = 

= j^[2Q[E{x),z][U{x,z)d{a)-V{x,z)e{-a)] - [U{x, z)S{a) + V{x, z)e{a)]] , 

(4.6) 

and it is obvious that, by construction, it obeys the homogeneous differential equation 

V'G'^(a;, ^) = 0. (4.7) 
Clearly, G^{x, z) has no support inside the future null cone of x. 
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Figure 4.1: Support of G"^''^'="\ G*^" and G 



Using the new Green's functions, the fields and can be defined. First, 
the singular field is 



ip^{x) = Airq J G^[x,z{t)] dr 



2(7 



and obeys Poisson's equation 



+ 



Tret 



2(7 



+1 

Tadv ^ 



''"adv 



V{x, z) (ir. 



(4.8) 



''"ret 



(4.9) 



It is noteworthy that the singular field does not depend on the entire past history 
of the moving particle, but only on its motion between the retarded and advanced 
proper times. 

The regular-remainder is defined as 



'q U{x, ;2)T^adv 



2(7 



- Q. 



J Tret 



Tret pTadv 

+ 9 , 

Tret 



V{x, z) dr, 



(4.10) 
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and by definition obeys the homogeneous differential equation 

VV^(x) = 0. (4.11) 

Since the field is a source-free field, it is smooth and differentiable at any point in 
spacetime and consequently at any point along the worldline of the particle as well. 
That is the most significant property of it, as far as the calculation of the self- force is 
concerned. 

To compare the scalar fields ^'^^'^ and ■0^, the difference of the two is calculated 
from Equations (4.10) and (4.2) 

^R(a;) _ ^^eif^^-) = / V{x, z) dr. (4.12) 

Using Equation (3.30) for the biscalar V{x, z) to expand the above integrand in the 
coincidence limit, the difference of the two fields becomes 

^ /"^adv 1 

^^(x) - r'^'^x) = -l [- -R{x) + 0{r)]dT 



2 Jrret 12 ^^^3^ 
= —qrR{x) + 0{r^), for x^F 

where the fact that the difference between the advanced and the retarded proper times 
is equal to 

Tadv - T,et = 2r + 0{r'), for X ^ r (4.14) 
was used. Taking the derivative gives 

VaV^ = VaV'^'^ + ^a\rR{x)\ + (tcrms that vanish as x ^ T). (4.15) 

The second term of the right-hand side of Equation (4.15) gives an outwardly pointing 
spatial unit vector and cancels the first term of the right side of Equation (4.3). Con- 
sequently, the fields '0^ and -0^®'^ give the same self-force. In addition, the averaging 
procedure that was necessary for the calculation of the self-force using i^f'^^'^ is not 
required when is used. 

A self-force calculation that has been demonstrated in [23] showed an additional 
benefit of calculating the self-force using the difference of the retarded and singular 
fields. The retarded field can be computed numerically and, since the singular field 
depends only on the motion of the particle between the retarded and advanced proper 
times, the knowledge of the entire past history of the particle is not required. That is 
a significant advantage of this method compared to the one that involves the direct 
and tail fields. 

Additionally, it is important that all the scalar fields used in this analysis are 
specific solutions of the homogeneous or inhomogeneous Equation (3.19). That makes 
them well-defined physical fields, a property that the direct and tail scalar fields do 
not have. 
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For those reasons, in this dissertation the self-force is calculated using 

= gV>^. (4.16) 

4.2 Electromagnetic Self-Force 

A detailed analysis for the self-force in terms of the direct and tail electromag- 
netic potentials was done in [14] and was presented in Chapter 3. In this section, the 
calculation of the self-force using the singular-source and the regular remainder fields 
is presented. 

It is assumed that a particle of electric charge q is moving on a worldline F : 
z'^ (r) in a background described by the metric gab- For the purposes of this section 
it can be assumed that there is no external electromagnetic field and the lowest order 
approximation to the particle's motion is a background geodesic. The particle creates 
an electromagnetic potential A" and an electromagnetic field F"-^, which interact with 
the particle, causing a self-force to act on it and making its worldline to deviate from 
a background geodesic to order q^. In the Lorentz gauge 

VaA- = 0, (4.17) 

the electromagnetic potential generated by the particle can be calculated by using 
Maxwell's equations 

and the different solutions of interest were described by the Green's functions of 
Section (3.2.2). 

If a potential A^^^^, analogous to ■0^^'^, is used to calculate the self-force in this 
case, its tail part is 

/Tret 
Vab'{x,z{T))u'' dr. (4.19) 
■oo 

Taking the derivative of it gives 

/Tret 
V,{Vab'u'')dT 
■oo 

1 1 X — Z 1 f^''"^ 

- gf {Rb'd' - ^9b'd'R) y" ^ - q / VcfKb'^i''] dr, 

z o r Ax^r 



(4.20) 



where Equation (3.53) for the bivector Vab' was used to derive the last expression. 
It is obvious that the same problems that were encountered in the scalar case are 
encountered here as well. The first term is not well-defined unless {Ra'b' " ^ga'b'R)u'^ 
is zero at the particle's worldline, the entire past history of the moving particle needs 
to be known and the potential used in the self-force calculation is not a solution of 
Maxwell's equations. 
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To overcome those difficulties, tfie singular Green's function in the neighborhood 
of the particle is defined analogously to the scalar case 



Gl,,{x,z)=G''^:;'{x,z) + ^Vaa'{x,z) 

= ^[Uaa'ix, Z)5{a) + Vaa'ix, z)e{a)]. 
OTT 



(4.21) 



It obeys the inhomogeneous differential equation 

V'Gl,{x, z) - R^'GUx, z) = -g-ha.'S'{x, z) (4.22) 
and gives the singular electromagnetic potential 

Al = ^7rJ {-9)-'^Gl{x,y)j\y)d'y, (4.23) 

which is a solution of the inhomogeneous Maxwell's Equations (4.18). By definition, 
both G^^/ and have no support inside the past and future null cone of x. The 
regular remainder is defined by 

A^{x) = ATix)-Al{x). (4.24) 

It has no support within the future null cone of x and it obeys the homogeneous 
Maxwell's equations. Consequently, it is smooth and differentiable everywhere in 
space, including every point along the worldline of the particle. The electromagnetic 
self-force can be calculated by the equation 

jr« = qg'^^VcAf - VbAf)z' (4.25) 

and is well-defined at the particle's location. It is stressed, again, that both A^ and 
A^ are well-defined solutions of Maxwell's equations. 

4.3 Gravitational Self-Force 

It is now assumed that a particle of mass m is moving on a worldline T : z" (r), 
in a background described by the metric gab- The particle causes a metric perturbation 
hab on the background metric. This perturbation obeys the harmonic gauge 

Vah""^ = (4.26) 
where hab is the trace-reversed version of the metric perturbation 

hab = hab - ^gabh%. (4.27) 
The linearized Einstein equations for it are 

V'hab + 2R:fh,a = -16777;;.. (4.28) 
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This metric perturbation interacts with the particle, causing a self-force to act on it 
and forcing its worldline to deviate from a background geodesic. 

If h^ab^, which is not a solution of Einstein's equations, is used to calculate the 
self-force, its tail part contains the bitensor Vabc'd' given in Equation (3.82). The 
contribution of this tail term to the derivative comes partly from the dependence of 
the proper time Tret on the point x. Specifically 



V7 self (tail) 
fKb --m 



— m 



Vahc'd'U'''u'^'V fTrA -ml V f\yabc'd'U'''u'^']dT 

ix—*T J —CO 



(4.29) 



Just like before, the first term is not well-defined at the particle's worldline, unless 
Rc'a'd'b'U'^ u'^ equals zero along the worldline. Also, knowledge of the entire past 
history of the particle is required in order to calculate the integral term. 

The singular Green's function Gabdd' {^i in the neighborhood of the particle 
can be defined as 

Glbc'd'{x, z) = GX"rf,(x, z) + -^Vabc'd'ix, z) 

^ (4.30) 

= -^[Uabc'd'ix, z)5{a) + Vabc'd' {x, z)0{a)]. 

OTT 

It gives the singular field h^f, which is an inhomogeneous solution of Equation (4.28) 
and has no support inside the past or future null cone of x 

hi = IGtt I (-g)-'^GlUx, y)T"'{y)d'y. (4.31) 

The regular remainder is defined as 

hl{x) = K^{x)-hl,{x). (4.32) 

It obeys the homogeneous differential equation 

V^h^, + 2R:fh^, = 0, (4.33) 

which means that it is smooth and differentiable everywhere in spacetime. Also, it 
has no support within the future null cone of x. 
Then the self-force can be calculated by 

^« = -m(^«^ + z''z'')z'z\VcKb - ^Vft/ifJ (4.34) 

which is well-defined at the location of the particle. It is important that both h^^ and 
are solutions of specific differential equations. 
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A significant conclusion can be drawn at this point. As befits the problem, the 

self-force is causing the particle to move on a geodesic of the metric {gab + ^afe)' ^'^ 
order m^. This metric {gab + h^b) itself, a homogeneous solution of the Einstein 
equations. An extended discussion of that point is given in [22] and [24] and is also 
presented in Chapter 8 of this dissertation. 



CHAPTER 5 

SINGULAR FIELD FOR SCHWARZSCHILD GEODESICS 



It was seen in Chapter 4 that the singular field is important in calculating the 
self-force. It is the part that must be subtracted from the retarded field to give the 
regular remainder, which is then differentiated to give the self-force. It is also clear 
that, since the singular field obeys the inhomogeneous Poisson, Maxwell or Einstein 
equations, it depends on the worldline of the moving particle and on the background 
spacetime. 

The coordinates that are used in this chapter to caclulate the singular field 
are the Thorne-Hartle-Zhang coordinates, abbreviated from this point on as THZ 
coordinates. Those coordinates were initially introduced by Thorne and Hartle [25] 
and later extended by Zhang [26]. A short discussion on them is presented in Section 
(5.1) of this chapter. 

The detailed calculation of the singular field of a scalar charge q for geodesies in 
a Schwarzschild background was presented in [23] . A brief discussion about that field 
is given here for completeness, since the regularization parameters associated with it 
are calculated in Chapter 6. The singular field is also calculated for geodesies in a 
Schwarzschild background, for a dipole generating a scalar field, for an electric charge, 
an electric dipole and a magnetic dipole each generating its own electromagnetic field 
and for a spinning particle generating a gravitational field. For all these calculations. 
Maple and Grtensor were used extensively. The Grtensor code is explained in 
Section (5.8) and presented in Appendix A. In Sections (5.3) through (5.7), only the 
results that the code gives are presented. 

5.1 Thorne-Hcirtle-Zhang Coordinates 

The Poisson, Maxwell and Einstein equations for the singular field assume a 
relatively simple form when they are written in a coordinate system in which the 
background spacetime looks as flat as possible. In the following, it is assumed that 
the particle is moving on a geodesic F in a vacuum background described by the 
metric Qab- Also, 7?. is a representative length scale of the background geometry, the 
smallest of the radius of curvature, the scale of inhomogeneities of the background 
and the time scale of curvature changes along F. 

A normal coordinate system can always be found so that, on the geodesic F, 
the metric and its first derivatives coincide with the Minkowski metric [18]. Such a 
normal coordinate system is not unique. The one used here is the THZ coordinate 
system and is used only locally, close to the worldhne of the particle. Specifically it 
is assumed that the background metric close to the worldline of the particle can be 
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written as 

9ab = Vab + Hab 



. (5.1) 



where rjab is the flat Minkowski metric in the THZ coordinates {t, x, y, z) and 
Also 



2Habdx"'dx'' = -£ijx'x^{df + 5kidx^dx^) + '^ekpqB'^y x' dtdx^ 



20 
21 



2 



5 



SijX X]^ p Sif^X dtdx -\- Xi€.jpq13 ^X^X ^pqiBj x^ 

o 



21 



1 



5' 



dx dx ^ 
(5.3) 



sHabdx'^dx' = 0{^)jx''dx' + 0{^)..dx'dx^. (5.4) 

In this and the following, a, &, c and (i denote spacetime indices. The indices 
i,j,k,l,n,p and q are spatial indices and, to the order up to which the calculations 
are performed, they are raised and lowered by the 3-dimensional flat space metric 6ij. 
The dot denotes differentiation with respect to the time t along the geodesic. Also, 
Eijk is the 3-dimensional flat space antisymmetric Levi-Civita tensor. 

If Hab consists of the terms given in Equation (5.3), the coordinates are second- 
order THZ coordinates and are well defined up to the addition of arbitrary functions 
of O(^). If Hab also includes the terms given in Equation (5.4), the coordinates 
are third-order THZ coordinates and are well defined up to the addition of arbitrary 
functions of O(^). 

The tensors S and B are spatial, symmetric and trace-free and their components 
are related to the Riemann tensor on the geodesic F by 

£ij — Rfitj (5-5) 
^ij ^ 2^i"^^pqjt- (5-6) 

They are of O(^) and their time derivatives are of 0(:^). 

For the calculation of the singular fields that follows, only the first two terms of 
Equation (5.3) are included in Qab- The remaining two terms are of 0{-^) and must 
be included in a higher-order calculation which must also take into account terms 
coming from the ^Hab part of the metric. 

The 'gothic' form of the metric is also defined as 



(5.7) 
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and its difference from the Minkowski metric is 

jjab ^ ^ab 

For the lowest order terms of Hab it can be found that 



where 



X-' 



2 10 

3 21 



^j^jX X"^ X ^ P I 

5 -I 



o 
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(5.8) 

(5.9) 

(5.10) 
(5.11) 

(5.12) 
(5.13) 



Using the simple symmetry properties of the tensors S and B, the following 
relationships can be shown for their components and the spatial THZ coordinates 

{x,y,z): 



and 



^ijk^S] + eikiB^j — €jkiB'l — 



( ^ijk^ rv^l ~l~ ^ijk^ln^ ^ikl^ n^^j ^jkl^ n"^ *^ 

{—eijkB\xi + eijkBinx'' — eikiB\xj + ejkiB\x^)x"'x^ 



= 
= 0. 



(5.14) 
(5.15) 



(5.16) 
(5.17) 



These relationships are used in the following sections when calculating the singular 
fields for different sources to simphfy the expressions for those fields. 

The calculation of the regularization parameters shown in Chapter 6 requires 
one to know the exact relationship between the THZ coordinates and the background 
Schwarzschild coordinates. A general procedure for finding that relationship for a 
given geodesic is to try to satisfy the three basic properties of the second order THZ 
coordinates: 

• On F: the coordinate t measures the proper time along the geodesic and the 
spatial coordinates x, y and z are equal to 0. Also, Qab = Vab and the first 
derivatives of the metric vanish there. 

• At linear, stationary order: H^^ = 0{x^). 

• The coordinates satisfy the harmonic gauge: dag"''' — 0{x'^). 
The details of this procedure are described in Appendix A of [23] . 
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For circular orbits in a Schwarzschild background of mass M, the THZ coor- 
dinates have been calculated in Appendix B of [23] using specific properties of the 
spherically symmetric background, which is much simpler than following the proce- 
dure just described. Their functional relationships with the Schwarzschild coordinates 
are given here in order to facilitate the discussion of the calculation of the singular 
fields and the regularization parameters. The Schwarzschild coordinates are {tg, r, 9, 0) 
and the Schwarzschild metric is 



ds^ = -(1 - '^)dtl + (1 - + r^de^ + sin^ 9dcj>^. 



(5.18) 



The circular orbit given by = ilts has orbital frequency equal to Q = (Mr~^)^, at 
Schwarzschild radius Tq. 

First, the functions X,y,Z, which are Lie derived by the Killing vector ^" = 



-M- + 0,^, are chosen 

ate 06 ' 



y 

z 



To sin(^) sin(0 — flt^ 
roCOs{9). 



(5.19) 

(5.20) 
(5.21) 



Then the functions x, y, z, also Lie derived by the Killing vector and the function 
t are defined in terms oi X,y,Z 



r sin 6* cos ( 
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1- 2M)^ 
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r2(i_ 2M)^ 
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[ - M^X'^ + y^{ro - 3M)(3ro - 8M) + 3Z^{ro - 2M)^] 



X^ 
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^-^A{rl-roM + ?,M^) 



8{ro - 2M) 



X^Z^ 



{2Srl - lUvoM + 123M2) + — — (Mr^ - 48roM + 33M2) 



14 



My^ 



4(r„ - 2M) 



(3^3 _ 74^2^ + 337roM2 - 430M^) 
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M'^y^Z'^{lro - 18Af) MZ^ 



56 



(5.22) 
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y = r sin ^ sin(0 — D,ts) ^ 



ro-2M\h 
- 3M 



My 



\r„-2M) 



Mxy 



(5.23) 



+ \2MX^iAro - 15M) 



+ :Vll4r^ - mroM + 89M^) + 2Z^(ro - 2M){7ro - 2AM)], 
2r;^(ro - 3m) 

+ [MX\l3ro - 19M) (5-24) 



To 

+ y^{Url - 36roM + 9M^) + Z^{ro - 2M){Uro - 15M)], 



3m,i rny 



^ To^ (1 - 

^ ^My 



'^Vo - + My + MZ^ 



r2^l-2M)^ro-3M)^ 2^° ' 3(ro - 2M) J (5.25) 

+ J^^(13r^ - 45r<,M + 31M2) + Z^(13ro - 5M)(ro - 2M)] . 
Finally 

X — xcos(Q'ts) — "UsinfO^to) 

+ + (5-26) 

y — xsm{frts) + ycos{frts) 



where fl^ = fljl - ^. 

There are two coordinate systems of interest. The first system is {t,x,y,z) 
which is a non-inertial coordinate system that co-rotates with the particle, meaning 
that the x axis always lines up the center of the black hole and the center of the 
particle. The y axis is always tangent to the spatially circular orbit and the z axis 
is always orthogonal to the orbital plane. As was already mentioned, for the spatial 
coordinates of this system it holds that: C^x = C^y = C^z = 0. 

The second system is {t, x, y, z) and is a locally inertial and non-rotating system 
in the vicinity of F. However, when viewed far away from F these coordinates appear 
to be rotating due to Thomas precession, as is clear from the terms involving the sine 
and the cosine of (fi^s) in Equations (5.26). It is also noted that: C^^p^ = C^z = 
>C^Vat = 0, but C^x, C^y and C^t are not equal to zero. This second system is used 
to calculate the singular fields. In order to avoid confusion of the field point x and 
the source point z with the THZ coordinates, the notation needs to be changed. The 
field point is denoted as p and its coordinates as x'^ and the source point is denoted 
as p' and its coordinates as x'^,. 
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5.2 Scalar Field of a Charged Pcirticle 

The singular field generated by a particle that carries a scalar charge q and 
is moving on a geodesic F : p'{t) (where, as usual, r is the proper time) of the 
Schwarzschild background must obey Poisson's equation 

V^V'^ = -47r^ (5.27) 
where the is written in the THZ coordinates and the source term is 

Q = q j{-9)-"^5\p-p\T))dT. (5.28) 
That singular field is derived in [23] and is equal to 

V'' = ^ + 0(^). (5.29) 

Here, instead of following the exact calculation of the singular field, which is performed 
in detail in [23] following a similar derivation in [27], 1 prove that it does indeed satisfy 
the scalar field equation to the order specified in Equation (5.29). 

The differential operator of the scalar field equation becomes, in THZ coordi- 
nates: 

If the field ip^ = (q/p) is substituted into this equation, the first term gives the 
expected 5-function singularity and the last two terms vanish since p does not depend 
on the time t. An explicit calculation shows that for the second term, x didjip is 
equal to zero. The remainder gives a term that scales as O(^). So 

V^V<^V„(^) = -ATrqS'ix^) + 0{^), for ^ ^ 0. (5.31) 

Consequently, for the remainder to be removed a term of 0{-^) must be added to 
q/p. So = g/p + 0(:^) is an inhomogeneous solution of the scalar wave equation 
and the error in this approximation is C^. 



5.3 Scalar Field of a Dipole 

The calculation of the singular scalar field generated by a dipole moving on 
a geodesic in a Schwarzschild background is presented in this section. The dipole 
moment is assumed to have a random orientation and its THZ components are denoted 
as Ka = {0,K^,Ky,K^). 

In this and the following sections, the subscripts (or superscripts) (0), (1) and 
(2) are used to indicate the order of significance of each term or component. The 
subscript (0) refers to the most dominant contribution, the subscript (1) refers to 
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the next most significant correction, which is calculated for the singular fields, and 
the subscript (2) refers to the next correction, the order of which is predicted for the 
singular fields. It is important to realize that multiplying two terms of order (1) does 
not necessarily give a term of order (2), because of the fact that there is no O(^) 
correction to the metric (see Equation (5.1)) and the fact that the first correction to 

2 3 

the metric is of O(^) while the second is of O(^). 

The scalar field of a dipole can be thought of as having the form 

*s ^ ^s^^ + ^s^^ + vi/S^) (5.32) 

and it obeys the differential equation 

V?o+i+2)^' = Vjo+i+2)(^(o) + + *(2)) = -47r^ (5.33) 
where the source term is given by 

g=-J K,V' [5\p - p'ir))] dr. (5.34) 

The zeroth-order term is the scalar field generated by a dipole that is stationary 
at the origin of a Cartesian coordinate system and it obeys the lowest order differential 
equation 

""(0)^(0) 

It is equal to 

.^s _ KjX' _ K^x + KyU + K^z 

(x^ + + ^2^3/2 

The first-order term obeys the differential equation derived from Equation (5.33) 



VL^fo) = -47r^- (5-35) 



V?o)*(i) + V?i)*?o)=0 Vfo)*fi) = -V?i)*fo) (5.37) 

which means that the V^j^-^ of the zeroth-order part of the field is the source term in 
the scalar differential equation for the first-order part of the field. In general, that 
source term is expected to contain the £'s and the B's and to give the first-order 
correction coming from the dipole's motion on the Schwarzschild geodesic. In this 
case. Equation (5.37) gives that 



' dz^ 

so the first order correction to the field can be set equal to zero. 

The source for the next order correction comes from the part of the that is 
of O(^) acting on the zeroth-order scalar field ^^{^y The differential equation is of 
the form 

d,d,%^ = O(^) X d,d^%^ ~ 0{^). (5.39) 
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That means that the next order term must be of O(^). 

Finally, the singular scalar field of a dipole moving on a Schwarzschild geodesic 
is equal to 

^,^ K^x + Ky + K,z ^^^p^ (5.40) 

5.4 Electromagnetic Potential of a Charged Particle 

In this section, the singular electromagnetic potential generated by a charge q 
moving on a Schwarzchild geodesic is calculated. Its general form is 

— ^s(o) + ^s(i) + ^s(2)- (5-41) 
Since in the vacuum background the Ricci tensor is 

Rab = 0, (5.42) 

this electromagnetic potential, as well as all those calculated in this chapter, obey the 
vacuum Maxwell's equations in curved spacetime 

^(0+1+2)^8(0+1+2) = -47rJ". (5.43) 
The zeroth-order term obeys the differential equation 

Vjo)Ag(o) = -47rJ" (5.44) 

the source term being 

r={qj {-9)"-S'\p - p'(r)] dr, 0, 0, o) . (5.45) 

The solution of this differential equation is the well-known Coulomb electromagnetic 
potential 

>ls(o) = (^, 0, 0, O). (5.46) 

The first-order correction to this electromagnetic potential obeys the differential 
equation derived from Equation (5.43) 

V?o)^s(i) = -V?i)A^(o). (5.47) 

The first-order part of (which contains the £^s and the S's) acting on the zeroth- 
order electromagnetic potential is the source term for the first-order correction. Sub- 
stituting the THZ components of ^s(i) 

^s(i) = (^s(i)>^s(i)>^s(i)'^s(i)) (5.48) 
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into the differential equation results in four differential equations, one for each one 
of these components. Each equation relates a specific sum of second derivatives of a 
component to a sum of terms of the form: 

qS X^X^X^X^ . , 

for the t-component, 

(^13 X' X' ^ ^ 

— '— — for the spatial components, 

where the dots denote appropriate indices. 

Solving these four equations is straightforward, once one notices that the solu- 
tion should have the form 

-SijX^x^ for the ^-component, 

q . (5-50) 

-e^^jB\x^x^ for the p-spatial component, 

each term multiplied by an appropriate algebraic factor. Substituting these expres- 
sions into the differential equations gives a set of simple algebraic equations for these 
factors, which can be easily solved to give the final expression for the first-order cor- 
rection: 

^S(l) = (~ 2p^*^^'^^' '2^^^^^^^^^^^^ ' Q^^^ij^^k^^^'^^ '2^^^ ij^^k^^ (5.51) 

The next order correction comes from the part of the that is of O(^) acting 
on the zeroth-order electomagnetic pothential. It gives a differential equation for each 
component of ^g(2) of the form 

a,a,.4^(2) = o(^) X a,a,Ag(o) ~ o(i^), (5.52) 

2 

which indicates that the next order correction must be of O(^). It is noteworthy that 
the first-order correction ^s(i) does not appear in the equation for the second-order 

2 

correction. That is because it only shows up in terms that involve the O(^) part of 
the metric, which are of the form 

O(^) X ~ 0(|j) (5.53) 

and must be included in a higher-order calculation. 

Finally, the singular electromagnetic potential of a charge q that is moving on 
a geodesic in a Schwarzschild background is equal to 

= f ([1 - lelB\x^x\ \e%,B\x^x') + O(^). (5.54) 
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5.5 Electromagnetic Potential of an Electric Dipole 

The calculation of the singular electromagnetic potential of an electric dipole is 
presented in this section. The dipole moment is assumed to point at some random 
direction and its THZ components are = {0,q^,q'^,q^). 

The singular electromagnetic potential can be written as 

= ^s(o) + ^s(i) + ^s(2) (5.55) 
and obeys the vacuum Maxwell's equations 

^(0+1+2)^3(0+1+2) = -47rJ" (5.56) 

where the source is 



J qiV'[S\p-p'{T))]V^dT, 0, 0, 0^. (5.57) 



The zeroth-order term is the electromagnetic potential generated by an electric 
dipole that is stationary at the origin of the Cartesian coordinate system, so it obeys 
the differential equation 

""■(0)^3(0) 

The solution to this equation is well-known and has only a t-component 



VLAg(o) = -47rJ^ (5.58) 



^s(o)=(^, 0, 0, 0). (5.59) 



The first-order correction to the electromagnetic potential obeys the differential 
equation 

^' (0)^3(1) = - V (i)^s(O)- 

If its THZ components are assumed to be 



^(0)^3(1) — ~'^(i)^s(o)- (5.60) 



^s(i) — (^s(i) > ^s(i) > ^s(i) ' "^s(i) ) (5.61) 

and are substituted into Equation (5.60), the differential equation for ^s(i) becomes a 
set of four second-order differential equations for these components. Each differential 
equation relates a sum of second derivatives of a component to a sum of terms of the 
form 

q£ X^X^X^X^X^ 

for the t- component, 

P (5.62) 
q 13 X' X' X' X' X' 

— — - — = for the spatial components, 

where, again, the dots denote appropriate indices. 
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These differential equations indicate tfiat tfie solution should be equal to a sum 
of the terms 

ft . ^ . 1 'T*^ 'yj rp^ i-Y^ . , rpj rp rp^ 



and 



(5.63) 



for the t-component and a sum of the terms 



X 



p6 p.i p.i 

e%q''B\x^x^x^ eijkq'BP^x'^XiX^ e^^^q^B{xPx''x^ 



(5.64) 



for the p-spatial component, each term multiplied by an appropriate mimcrical co- 
efficient. In fact, using Equations (5.15) and (5.17), the last two terms that are 
expected to show up in the solution for the p-component can be eliminated in favor 
of the remaining four. Substituting these expressions into the differential equations 
gives simple abgebraic equations for the coefficients. The final expressions for the 
components of ^s(i) 



— — 
^S(l) - o 



1 qiSjkX^x^x'' 



^s(i) - 



e^^^q'^B\x''XkX^ e^^^q^B\x^ xix^ e-^^jq'^B\x^ x^x^ 



(5.65) 



2p3 



+ 



2p3 



+ 



2p3 



The order of the next term in the expansion of the electromagnetic potential 
can be predicted. It is the solution of the differential equation whose source term 
comes from the 0{^) part of the acting on Ag^Qy That differential equation has 
the form 



Consequently, the second-order correction must be of 0{:^). The terms that involve 
the first-order correction ^s(i) '^'^ ^'^^ contribute to the equation for the second-order 

correction, because they involve the O(^) part of the metric and that results in terms 
ofO(^). 

Finally, the singular electromagnetic potential for an electric dipole moving on 
a Schwarzschild geodesic is equal to 



/I* - — 



4 



1 



e%q'B\x'xkx'' + e^q'^B'kX' 



^^x^xix' + e%q'B\x^x'x']+0{:^^: 



(5.67) 
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5.6 Electromagnetic Potential of a Magnetic Dipole 

In this section, the singular electromagnetic potential of a magnetic dipole mov- 
ing on a geodesic of the Schwarzschild background is calculated. The magnetiza- 
tion m" is assumed to point at some random direction and its THZ components are 
m" = (0, m^, m^, m^). 

The singular electromagnetic potential can be written as 

= ^s(o) + ^s(i) + ^s(2) (5.68) 
and obeys the vacuum Einstein equations 

V?o+i+2)^s(o+i+2) = -47rJ" (5.69) 

where the source term is 

J° = 0, 

j9^y" ^qijQ. mj8^{p-p\r)) yf^dT. 



(5.70) 



The zeroth-ordcr term is the electromagnetic potential generated by a magnetic 
dipole that is stationary at the origin of a Cartesian coordinate system and is the 
solution of the differential equation 

V?o)^s(o) = -47rJ". (5.71) 

Its THZ components are 

^S(O) - ij—^^^ ij—^^^ i3~^)- V^-'^J 



The first-order term obeys the differential equation derived from Maxwell's 
Equations (5.69) 

^■(0)^8(1) = "'^(1)^8(0) 

which indicates that the source term for ^g(i) is the V^^-j of the zeroth-order part of 
the potential. The THZ components of ^s(i) 



^(01^8(1) — ""^(1)^8(0) (5.73) 



^s(i) — (^s(i)'^s(i)'^s(i)'^s(i)) (5-74) 

and when substituted into Equation (5.73) the result is a set of four second-order 
differential equations, one for each one of those four components. Each equation 
relates a sum of second derivatives of one component to the source term which consists 
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of terms of the form: 

m B X X X- 



^5 



for the t-component, 



TT.'r.'r: 

for the spatial components, 



c (5.75) 
m c XXX X X' 



where the dots denote the appropriate indices. 

Solving these differential equations is tedious but not difRcTilt. because each 
equation involves only one component of ^s(i) and only one of the two tensors £ and 
B. A careful look at the equations indicates that the solution should be a sum of the 
terms 

rUiB'^x^XkX^ rriiBjkx'x^x^ ,^ 
^— and (5.76) 

for the ^-component and a sum of the terms 

e^jTYi^ £\x'^ xix^ e^^^m^Sux-'x'^x^ e^^jUi'^S^x^xix'' 



p3 ' p3 ' p3 

^^jm^£\x^XkX^ eij^m^'EP^x^xix^ eijk'm^S^x^x^x'P 
p3 ' p3 ' p3 



(5.77) 



for the p-spatial component, with appropriate numerical factors in front of each term 
so that the equations are satisfied. Using Equations (5.15) and (5.17), the first and 
last terms expected to show up in the final expression for the p-component can be 
eliminated, since they can be expressed as linear combinations of the remaining four 
terms. Substituting these expressions into the differential equations gives a system of 
four algebraic equations for those factors. Solving these algebraic equations is trivial. 
The result is that the first-order components of the electromagnetic potential are: 

^S(l) ^ ~3 \^7:1TT'i^jkX^X-' X — —m^BijX^X Xk\ 

\ . . 1 1 (5-78) 

^s(i) = Y^jVPi-Eux^x^x^ - -^^■mkE\x^x^x^ - -eijkm'EP^x''xix^] . 

P Zi Zi 

The order of the next correction to the singular electromagnetic potential can 
be predicted. It is the solution of the differential equation that has the O(^) part 
of the acting on the zeroth-order electromagnetic potential as the source term. 
Specifically, it looks hke 

a,a,.4gp^ = o(^) X a,a,^g(o) - o{^\ (5.79) 

meaning that the ^s(2) correction is of O(^). As in the previously studied cases, the 
terms that involve the first-order correction do not contribute to this equation, 

since they involve the O(^) part of the metric which results in terms of O(^). 
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Finally, the singular electromagnetic potential for a magnetic dipole moving on 
a Schwarzschild geodesic is equal to 

XX 2 
Al ^ ^[-rriiBjkx'x'x^ - -m'B^,x' x^Xk] + O(^) 

p p m^x^ r p m^Skix^x'^x^ 1 p rnkS^iX^x'^x^ 1 rrfS^^x^xixK p 

(5.80) 



5.7 Gravitational Field of a Spinning Particle 

The calculation of the singular gravitational field of a spinning particle mov- 
ing on a Schwarzschild geodesic is presented in this section. The particle is as- 
sumed to have a small angular momentum pointing at some random direction A"- = 
(0, A^, Ay, A') in THZ coordinates. 

The singular gravitational field can be written as 

and obeys the linearized Einstein equations 

V(o+i+2)4°^6^^^^ + '^^{0+1+2) ab^^^f^d^^^ = -IfiTrTab (5.82) 
where ^ 

ab — hsab— ^dab^S c (5.83) 

is the trace-reversed version of hsab- 

The zeroth-order part of the singular field is the gravitational field generated by 
a particle with angular momentum A" that is stationary at the origin of a Cartesian 
coordinate system. 

For the angular momentum pointing along the z-axis, that is the well-known 
Kerr solution with the mass set equal to zero. Since the angular momentum is assumed 
to be small and the effects of the mass of the particle are not taken into account, 
only the terms of the Kerr metric that are linear in the angular momentum need to 
be considered. In the Boyer-Lindquist coordinates {t^h, f^, <P) around the spinning 
particle, h^^b (where the superscript Z denotes that this is the part of the zeroth-order 
gravitational field coming only from the ^-component of the angular momentum) is 
equal to 



hl,,^--^A^si^'e 



/ 1 \ 





\ 1 / 



(5.84) 
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Using Grtensor, this expression can be easily converted into the equivalent expres- 
sion in THZ coordinates 



f^S ab - ^3 



/ 





y 


—X 


0\ 




y 













—X 











V 














(5.85) 



The relationships between the Boyer-Lindquist and the THZ coordinates used for 
this conversion are the usual relationships between the spherical and the Cartesian 
coordinates, namely 

t — tBh 

r — y^x^ + + 

y/^^+F (5-86) 

= arctan 

z 

(f) — arctan — . 

X 

This is sufficient because the corrections to these relationships that involve the angular 
momentum would give terms of higher order in the angular momentum and must be 
ignored in this analysis, since only first-order terms in the angular momentum are 
kept. 

The components of the angular momentum along the x and y axes must be 
treated separately, because of the axial symmetry of the Kerr metric. The analyses 
for the angular momentum being along the x and y axes are very similar to that 
for the angular momentum being along the z axis and the only change comes from 
the different form that the relationships (5.86) have. Specifically, when the angular 
momentum points along the x axis, the axial symmetry is around the x axis and the 
relationships used are 

t — tBh 

r — y^x^~-\-y^~+~z^ 
, = arctan 

X 

(p — arctan -. 

y 
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For the angular momentum pointing along the y axis, the axial symmetry is around 
the y axis and the relationships are 



t = t 



BL 



9 — arctan 



(5.88) 



arctan ■ 



X 



Adding all the contributions that result from this analysis, the zeroth-order singular 
gravitational field becomes 



"-S ab 



( 

{-A^z^A^y) 

{A^'z - A'x) 
\ {-A^y + Ayx) 



{-Ayz + A'y) {A'z - A'x) 









-A'^y + Ayx) \ 



/ 



(5.89) 

The first-order correction to this gravitational field obeys the differential equa- 
tion derived from Equation (5.82): 



Y72 I OF? cdi,(l) 

^ (0)"S ab + ^-^(0) ab'^Scd 



^(l)"Sa6 ^-"(l)a&"S 



\cd 



(5.90) 



SO the first-order and the first-order Riemann tensor acting on the zeroth-order 
solution give the source for the first-order correction. The first-order correction must 
be a symmetric tensor so it is assumed to be equal to 



"S ab — 



'Hx 



'Hx 


Iky 




\ 


ll-XX 


Ii'xy 


ll'XZ 




ii'xy 


Ibyy 


Ibyz 




ll'XZ 


Ibyz 


ll-ZZ 


/ 



(5.91) 



Substituting it into the first-order equation results in 10 differential equations. There 
is one set of four differential equations for the four diagonal components, each equation 
containing all four diagonal components. There is also one differential equation for 
each one of the t — i components and one differential equation for each one of the i—j 
components, for i ^ j. In each equation, a sum of second derivatives of components 
is related to a sum of terms of the form 



A B X^X^X^X^X^ 



A £ x x x-x'x- 



, for the t — t and p — q components. 



, for the t — p components. 



(5.92) 



where the dots denote the appropriate indices for each term. 
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Solving the differential equations in this case is slightly more complicated than 

in the previous cases, mainly because of the fact that four of them involve all diagonal 
components rather than only one of them. Still, the process becomes significantly 
easier if one notices that the t — t component must be an appropriate sum of terms 
of the form 

and —„ , (5.93) 

each t — p component must be a sum of terms of the form 



p3 ' p3 ' p3 

,Ai£',.r'x'.r'' f,jiA'&'/x,x'' t,nA'&'v>>:j^x'' 



(5.94) 



pi, pi, pi, 

and each p — q spatial component must be a sum of terms of the form 

^p*^qk^ ^q*^pk'^ -^l^ ^p*-^kl^q^ ^ ^q*^kl^p^ ^ 

p"^ p-^ p^ p"^ 

■^k^pq^ -A JSf^iXpXqX Tjpq-A.i]3ji^X X"^ X TjpqA ^ijX^ X}^X 

P*^ P-^ P"* P*^ 



(5.95) 



A^BpkXqXlX^ A^Bq};XpXlX Aj;BplXqX X A^BqlXpX 



k^l 



p3 ' p3 ' p3 ' p3 



with appropriate numerical coefficients in front of each term. Equations (5.15) and 
(5.17) can again be used to eliminate the first and last terms in favor of the remaining 
four, for the expression for the t — p components. Substituting these sums into the 
differential equations gives algebraic equations for the coefficients, which are fairly 
easy to solve. The result is that the first-order correction to the gravitational field 
has components 



''Stt — 


J^-i^jj^X X^ X 






1 

7 






"stp — 


^P'ij'^'l^ k^ ^ ^ ~ 


^pij-^^lk^"^ ^ ^ ~l~ ^^ijl-^^J^ ^k^ 




1 

7 


^.^.'i^SpqX X jX ~|~ 


2 -2 
~-A ^ijXpXqX"^ ~\~ ^i(^pXq'^XjX'^ ~ 



10 ■ ■ 2 ■ ■ 2 ■ ■ k ' k 

"^-^{p^qji-^ •^j-^'^ '^A(j,Xq^l3ijX X'^ -\- —Tjpq(^Ajl3jf;X X'^X A BjjX"^ X f^X ^ 

(5.96) 

The order of the next correction can be predicted. That correction is the solution 
of the differential equation 

^W^slb + 2-R(0) ab^sld = - V(2)4°a6 ~ '^^(2) a b^sld (5.97) 
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where V(2) and R^2)ab come from the O(^) part of the metric. So the equations 
look hke 

^3 1 

(5.98) 



and the second correction to the gravitational field is of 0(:^). Again, h^^^j^ does not 

2 

show up in the equations for the second correction, since it relates to the O(^) part 
of the metric and results in terms of O(^). 

Finally, the singular gravitational field due to a spinning particle moving on a 
Schwarzschild geodesic is equal to 



hstt ■ 
hstp 



h: 



Spq 



-2 



+ 



-2 



. . 2 • • • -2 j • 

10 .2 -2 . . ^ . . ^ 

~^^[p^q)iX XjX"^ —jA.(j)Xq^^ijX X'^ ~\~ —7]pq(^Ai^jj^X X'^X A ISijX'^ Xj^X ^ 



(5.99) 



5.8 Grtensor Code 

The Grtensor code (running under Maple) used to derive the differential 
equations for the first order correction to the singular field is given in Appendix A. 
Since the case of the gravitational field is the most complicated one, the analysis is 
presented for the gravitational field generated by a spinning particle with the angular 
momentum pointing along the THZ z-axis. The analyses for the scalar fields and the 
electromagnetic potentials are very similar, and can be easily deduced from that for 
the gravitational field. 

An effort was made to keep the symbols in the code in accordance with the ones 
used in this chapter for the various quantities. In the situations where that is not the 
case, the comments in the code should make the notation clear enough for the reader 
to follow. 

The parameter e is used to keep track of the order of each term in the com- 
ponents of the tensors S and B and is set equal to 1 at the end. Throughout the 
calculation, only first order terms are kept. Specifically, the Christoffel symbols and 
the components of the Riemann and Ricci tensors are calculated first and all their 
terms that are of order higher than 1 are set equal to zero. Doing that makes the sub- 
sequent analysis significantly simpler and the running time of the code significantly 
shorter. 
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After the various quantities associated with the problem are calculated, the 
test tensor hbartest (a,b) , whose exact dependence on the THZ coordinates is not 
specified, is used as a trial solution in the linearized Einstein equations and the source 
term coming from the zeroth-order solution hO(a,b) is examined. That helps identify 
the terms that should be expected to show up in each component of the solution 
hbar(a,b). Specifically, it helps determine which tensor's components, £ (denoted as 
EE in the code) or B (denoted as BB), should show up in each component of hbar(a,b) 
and gives an idea of how they should be contracted to the spatial THZ coordinates 
X, y, z. The terms that result from this analysis are multiplied by algebraic factors 
and the appropriate sum is substituted into the linearized Einstein equations. The 
result is a simple system of algebraic equations which can be easily solved to give the 
values of the algebraic factors. That completes the solution. 

The last component of this analysis is a simple confirmation performed for the 
gravitational field that was calculated, that the algebraic coefficients obtained do, 
indeed, give the required solution. The confirmation is simply done by replacing the 
initially unknown algebraic coefficients with their exact values in the expression for 
the solution and substituting that expression into the Einstein equations. Despite the 
fact that it was not explicitly mentioned in Sections (5.3)-(5.7), that confirmation was 
performed for all singular fields and potentials that were calculated. 



CHAPTER 6 

REGULARIZATION PARAMETERS FOR THE SCALAR FIELD 



A mode-sum regularizarion procedure for the scalar singular field is presented 
and implemented in this chapter. This regularization procedure was first proposed by 
Barack and Ori in [28] , where they described the calculation of the regularization pa- 
rameters for the direct part of the self-force on a particle carrying a scalar charge. The 
procedure was later implemented by different groups for the calculation of the regular- 
ization parameters for the direct part of the self-force on a scalar charge on different 
geodesies [28-33] and also for non-geodesic motion [31, 33] around a Schwarzschild 
black hole. The calculation of the regularization parameters has also been performed 
for the direct part of the electromagnetic self-force [34] and for the direct part of the 
gravitational self- force [34-36], for arbitrary geodesies around a Schwarzschild black 
hole. 

Even though the regularization procedure was initially described for the contri- 
bution of the direct part of the scalar field to the self-force, it can be used equally 
successfully for the contribution of the singular scalar field to the self-force, as was 
demonstrated in [23]. In that paper, the regularization parameters for the self-force 
on a scalar charge in circular orbit around a Schwarzschild black hole were calculated 
and the self-force results ended up being in excellent agreement with the results that 
were derived using the direct scalar field [28-31]. 

This chapter begins with an outline of the regularization procedure for the scalar 
self-force. The description closely follows that given by Barack and Ori [28] for the 
direct self-force but is presented here for the singular self-force instead. Then, the 
regularization parameters are calculated for the singular scalar field (rather than the 
scalar self-force) of a charged particle that moves on an equatorial circular orbit in 
a Schwarzschild background and the results of [23] for the scalar self-force are repro- 
duced. Finally, the regularization parameters for the first derivative of the singular 
part of the self-force are also calculated. 

6.1 Regularization Procedure 

As was shown in Chapter 4, the self- force on a particle that carries a scalar 
charge q can be calculated from the equation 

JT^ = q lim V„^^ = q lim VaU'""^ - (6.1) 

where p' is the point along the worldhne of the charged particle on which the self-force 
needs to be calculated and p is a point in the vicinity of p'. 
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It is assumed that the charge q is moving in a Schwarzschild background of 

mass M and the Schwarzschild coordinates are {ts,r,6,(f)). For the calculation of the 
retarded field, the source term in Poisson's Equation (3.19) can be decomposed in 
terms of spherical harmonics and the retarded field can be written as 

oo I 

ij'-'\t,,r,e,<P) = J2Y. i^Zir,QYi^{e,<P). (6.2) 

1=0 m=-l 

Then, the /m-components of ip^'^^ can be calculated numerically. That calculation 
is discussed in great detail in Chapter 7. Here it is just noted that the important 
property of the V'fm'^ their first r-derivatives is that they are finite at the 

location of the particle, even though is singular there. 

If the spherical harmonic decomposition of the singular field is also considered, 
Equation (6.1) becomes 

oo I 

= q hm, W - ^f-) = ^ li"^, E X - O Yim. (6.3) 

1=0 m=—l l,m 

It is helpful to define the multipole /-modes of the two contributions to the self-force, 
which result after performing the m-summation of each term individually in Equation 
(6.3). Specifically 

m 

which gives for the self-force 




In Equation (6.5), the difference in the multipole /-modes must be taken before the 
summation over / is performed. 

From this point on, the discussion of the regularization procedure becomes spe- 
cific to the problem of the scalar field (q/p), since more detailed results are available 
for this case. However, a similar analysis can be done for any other scalar field. The 
goal is to find a function hia such that the series 

J2(^a'-hla) (6.6) 

I 

converges. When such a function is found, the self-force can be written as 



I 



(6.7) 
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where 

E,= \iv^^{Ti-hia). (6.8) 

I 

Because of its definition, the function hia should be calculated by investigating the 
asymptotic expansion of for large On the other hand, because the self- force is 
known to be well-defined, !F12^ and Tf^ are expected to have the same large-Z behavior, 
so hia can be determined by the asymptotic behavior of JF^ instead. 

The singular part of the self-force consists of terms of different order in the 
limit p — > p' and it has been shown [23] that, in principle, only the first three of those 
terms are expected to give non-zero contributions, for the field {q/p)- However, for 
reasons that will become clear shortly, the next order terms are included and 

The superscripts A,B,C,D indicate the different orders, A coming from the most 
dominant term, B from the next more dominant and so on. The superscript E 
refers to all terms of order higher than that of J-'a^^\ The mode-sum regularization 
procedure amounts to performing the spherical harmonic decomposition of each such 
term, which results in an expression of the form 

= E(^a" + + + + ^^)^- (6.10) 

l,m 

with A^^ corresponding to J^^'^\ etc. For the simple case of a scalar charge moving 
in a Schwarzschild background and generating the scalar field {q/p), the parameters 
-B^™, C^™ and D]^ have been shown [23] to have a very simple form so that, when 
the explicit expression for the spherical harmonics Y/m is substituted into Equation 
(6.10), the summation over ni can be performed and the result is an expression of the 
form 

^ - E [Ml + i) + i^. + C„ + (2,-i)p, + 3) + 

where the regularization parameters Aa, Ba, Ca, Da and Ea arc /-independent quan- 
tities which do depend on the background geometry and the characteristics of the 
orbit. 

Finally, the self-force can be calculated by 

oo 

' »— >»' 

'=0 (6.12) 

One important point that should be made is that the infinite sum over I must be 
performed, in order for the self-force to be calculated. Notice, however, that the 
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contributions for large I get less significant for the terms containing Da and Ea. The 
reason for including these terms is now clear. Even though the sum over / of each of 
these two terms is exactly equal to 0, including these terms improves the convergence 
of the sum. An additional benefit of including these terms is that the approximation 
to becomes more differ entiable, as is explained in [23] . 

6.2 Order Calculation of the Scalar Field 

It should be obvious from the analysis of Section (6.1) that, in order to calculate 
the regularization parameters for the singular field generated by a charge q, it is 
necessary to have an expression of (1/p) in which the order of each term is known. 
The derivation of such an expression is presented in this section for an equatorial 
circular orbit of radius To in the Schwarzschild background. It is noted that the 
results were derived using Maple extensively. It is also noted that the derivation was 
presented in [23] where the results of it for the radial derivative ^r(~) were given. For 
simplicity, the scalar charge q is set to 1. 

In order for the calculation of the self-force regularization parameters to be made 
easier, the Schwarzschild coordinates can be rotated, as explained in [30]. Specifically, 
new angles © and $ can be defined in terms of the usual Schwarzschild angles by the 
equations 

sin ^ cos(0 — fits) = cos© (6.13) 
sin 6* sin(0 - Q^s) = sin© cos $ (6-14) 
cos^ = sin© sin$ (6.15) 

so that the coordinate location of the particle is moved from the equatorial plane, 
where 6* = |, to a location where sin© = 0, for a specific tg. Such a coordinate 
rotation preserves the index I of any spherical harmonic Yim{9, 0). That means that 
any Yim{0, (t>) is mapped into a linear combination of spherical harmonics Yim'{Q-i $), 
where m' = —I, Consequently, each /-multipole mode of the field or the self- 
force that results after summation over m is the same, regardless of which angles, 
{9, 4>) or (©, $), are used for calculating it. 

The benefit of this coordinate rotation for the calculation of the regularization 
parameters for the singular part of the self-force can be understood if one remembers 
that in the limit p — > p' the angle © is equal to 0. That means that Yijn{0, $) has to 
be used, for which 

Yi^{0,^) = 0, form^O (6.16) 
y,„(0,$) = ^J?±^, iorm = 0. (6.17) 

So the sum over m can, after this coordinate transformation, be replaced with just 
the m — term. However, for the regularization parameters of the singular field, 
the limit p ^ p' must not be taken, so the Yi^^s for all m's must be taken into 
consideration. The m — spherical harmonic must be considered only when the 
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regularization parameters for the self-force are derived from the ones for the scalar 
field and when the regularization parameters for the first r-derivative of the self-force 
are calculated. 

A comment on the order of each term needs to be made at this point. The 
parameter e" is used to indicate a term of order in the coincidence limit p — > p'. 
In that limit, r — > To and © — > 0. That means that the factor (r — Tq) is of order e 
and the factor (1 — cos 0) is of order e^. At the end of the calculation, the parameter 
e can be set equal to 1. 

The relationships between the Schwarzschild coordinates {tg, r, 9, (p) and the THZ 
coordinates {t, x, y, z) for a circular orbit on the equatorial plane of a Schwarzschild 
background that were given in Section (5.1) are used for this calculation. As was 
mentioned earlier, = + y'^ + in terms of the spatial THZ coordinates x, y and 
z. However, it is clear from Equation (5.26) that the relationship 

x"^ + y"^ ^ x'^ + f (6.18) 

holds between {x,y} and {.t, y}, so the sum (5^ + y'^) is used to calculate p^. The 
Equations (5.22), (5.23) and (5.24) are substituted into the expression for p^, with the 
Schwarzschild angles 9 and replaced by the new angles © and That substitution 
gives that the lowest order term for p^, denoted by p^, is 

p2 = — —rT7 + 2rl- — x(l-cos©, 6.19 

where 

A = r - To (6.20) 

and 

X^l — sin2$. 6.21 

Clearly, p^ is of order in the coincidence limit, as should have been expected. Next 
the variables 0, $ and r are eliminated from the order expression of p^ in favor of 
the variables A, p and by using Equations (6.19), (6.20) and (6.21). Finally, the 
result is inverted and the square root is taken, in order to obtain the order expression 
of(l/p). 

The result of this calculation is a very long expression. Here, I only give the 
terms of this expression that are necessary to calculate the regularization parame- 
ters of the singular field, the singular self-force along the radial direction and the 
r-derivative of the singular self-force along the radial direction. Exactly how that is 
determined will become clear shortly, when the general term A^p^x"^, where a, p 
and q stand for integers, will be discussed. 
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The terms of interest are 



P P 



To-SM 1 _ 1 

2ro(ro -2M)x~~o 



A 

— + 
P 



2ro - 3M 



ro-3M 1 



To - 3M 
8r2(ro - 2M) 



1 To+M 1 



2ro(ro - 2M) 2r<,(r<, - 2M) x, 
P 



+ 



2ro - 3M brl - 22roM + 21M^ 1 Sr^ - 22roM + 21M2 1 



+ 



2r2(r, - 2M) 

M{ro - 2M) {to - M) {to - AM) 1 
2rt{ro - 3M) 8(ro - 2M)rt X 



P 



+ 



^,A3 A^ AM 
^ p p"^ p' ) 

+ 0{e% 



16r^{ro - 2Mf x^ 



Ap 



(6.22) 



6.3 Scalar Monopole Field 

As is clear from the previous calculation of (1/p), the angular dependence of 
every term shows up in factors of the form f^x~'^) where p is an odd integer and q = 
0, 1,2,.. .. The spherical harmonic decomposition of that factor, which is necessary 
when calculating the regularization parameters for the scalar field, is given in detail 
in Appendix B. The result is 



P'X 



2rl{r, - 2M) 



3M 



p/2 oo I 
1=0 m=-l 



l.m 



M 



2M 



where 



7 



3M) 



(6.23) 



(6.24) 



2ro{ro - 2M)2 

and the coefficients Ef'^{'y'^, ^ ^2m ) given in Equation (B.41). It is stressed that 



,2s 



'-)"" and the term (7^) 2 



all the r-dependence of p^x resides in the sum Y1T=q 
Equation (B.41) for ^'f^, and is always proportional to powers of A = (r 



hn+l 



m 
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At this point, a note on the term (7^) 2 for an odd integer /3 is in order. For 
that term, the square root of 7^ must be considered. One might think that 



3M 



2ro{ro - 2M)2 



3M 



2ro(r„ - 2M)2 



(6.25) 



But that imphes that being on the equatorial plane with = and approaching the 
particle by taking the limit r — > To, would give for the leading term p in the expansion 
of p: 



2r 



2 ' o 



2M 



Wo-3M 



2r 



,ro-2M 



ro-2M 



X2 ( — + 1 — cosO 



(7^ 



A. 



X 



(6.26) 



According to this, the sign of the leading order term of p could be either positive or 
negative, depending on whether r > Tq or r < Tq. That is clearly not correct, since by 
definition p = ^/ x'^ + y'^ + z"^, which is always positive. For that reason, taking the 
square root of 7^ always implies 



To - 3M 
2ro{ro - 2Mf 



(6.27) 



For the self-force along the radial direction to be calculated, the r-derivative of 
(1/p) and the limit r — > Tq have to be taken. That means that any term of order 
(r — or higher gives, after the limit is taken, no contribution to the self- force. 
However, as has already been mentioned, it is desired to calculate the first derivative 
with respect to r of the self-force, which is the second derivative with respect to r 
of the field, in the limit r — > Tq. Consequently, terms of order (r — have to be 
retained, because they do give a contribution at that limit, while any terms of order 
(r — To) 2 or higher can be disregarded. 

As is clear from Equation (B.41) of Appendix B, the general term N^pp'x^^ in 
Equation (6.22) has two pieces that contain the r-dependence. The first piece comes 
from the sum over s in and gives terms proportional to 7^+2* . Such a term should 
be kept only for a -|- 2s < | or 

(6.28) 



5 



a 
2' 



Since 

a > 



s > 0, the 7°'+2s contribution can be immediately disregarded for terms with 
_ |, while for terms with a < | a limited number of values of s have to be 
retained in the sum. The second piece that contains r-dependence comes from the 
term (72)f+"+i and is proportional to ^p+"+2n+2 Such a term should be kept only 
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for p + a + 2n + 2 < I or 

n < - - (6.29) 

Since n > 0, all terms of Equation (6.22) for which p + a > | do not contribute 
through the ^p+"+2n'+2 piece. Combinations of p and a for which p + a < ^ should be 
examined individually and a specific number of n's must be retained. The terms of the 
expansion (6.22) for (1/p) for which just the order (and not the explicit expression) 
is given are terms that fall in both categories that according to this analysis can be 
ignored. 

I present now the calculation of the regularization parameters coming from the 
terms of different order, for the singular field (1/p). The abbreviation used for the 
hypergeometric function is 

1 TYl A/f 

F.^.™..F.(a,A+-;U + i;__), (6.30) 

and it is also noted that the hypergeometric function 

F°'° = 2i^i(a,-;1; ^.J (6.31) 

is denoted by Fa, as is done in [23]. Also, for the three sums over k, v and A the 
abbreviation 



E 



2Z + 1 r(/-|m| + l) 

TT [r(/-m + l)r(Z + m + l)]^ 



■X 



^ V{k + -2k- \m\ + 1) ^ r{u + l)r(M -jy+i) 

r(^) 



r(A + i)r(^-A + i) 



A=0 

(6.32) 

that is used in appendix B is also used here. In addition, the upper limit of the 
n-summation is denoted by 

N^l-\m\-2k + 2v. (6.33) 
For the lowest order term of Equation (6.22) 

T(_i) = e-'^ (6.34) 

the exponents of A"p^x~^ are a — 0,p— —1, q — 0. Equations (6.28) and (6.29) give 
that only the s — and s = 1 terms of the first piece and only the n = term of the 
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second piece of Ei ^'^ need to be kept. Consequently 



l,m 



(6.35) 



where, using Equation (6.23) 

To - 3M V 
2r2(ro-2M)_ 

Tr, - 3M 



E, 



-1,0 



2rl{ro - 2M) 



A.m 



(6.36) 



iV 



n=0 



r_^y+n22n+i r(A^ + 1) 

(2n + l)!! r(iV-n + l) 



2F, 



A,m 
1/2 



- 2)^3/2 7 



or, when the exphcit expression for 7^ and its square root are substituted 



A 



Im 



3M 



3M 



2r„(r„ - 2M)2 



l^l^i 



n=0 



_^y+n22n+i r(iV+l) 

(2n+l)!! r(iV-n + l) 



3M 



V2 2^ 3/2 2r„(r„-2M)2' 



(6.37) 



As is explained in Appendix B, only the even m's should be included in the sum of 
Equation (6.35). 

The regularization parameter given in [23] for the self-force can be easily 
obtained from this result. First the m = parameter is considered, for reasons that 
were explained earlier. That makes the sums over A and over v equivalent to the 
terms with A = and v — For the hypergeometric function F^'^ it can easily be 
proven that 



To - 2M 



(6.38) 



ro-SM 

using Equations 15.3.3 of [37]. Then the first derivative of Ai^ with respect to r is 
taken. That makes the term proportional to F^^^ in Equation (6.37) vanish, while the 

term proportional to F^j^ gives a factor of A which vanishes when the coincidence 
limit is taken. For the first term the derivative of |A| gives sgn(A). Finally, the 
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spherical harmonic Yio{0, is substituted with \/%^- This procedure gives 



d_ 
dr 
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1=0 
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1=0 
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E(-i)'2 
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k=0 



T{k + l)T{l-2k + l) 



(6.39) 



The sum over k can be easily calculated for any value of I, using Maple. A general 
proof that it is equal to ^/^^ for any / cannot be given. However, for every value of I 
that was tried, the sum ended up being s/tt, which gives 



-sgn(A)) 



p^p' 



[rojro - 3M)] 
r2(r„ - 2M) 



1 oo 



(6.40) 



1=0 



which is the result of [23]. 

Now, the second derivative of the term T(_i) with respect to r is taken, which 
is equal to the contribution of that term to the first derivative of the self-force. The 
m — component is considered since the coincidence limit must be taken, so A = 
and 1/ — as well. Only the term proportional to F^j^ survives the differentiation, 
since it is proportional to A^. The result is 
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(6.41) 



The zeroth-order term contribution that is considered in Equation (6.22) is 



7^(0) = e° 
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. (6.42) 



For the piece that is proportional to Ap~^ the exponents are a = l,p = —1, = 1 for 
the term that contains x ^i-nd a = 1, p = —1, g = for the term that docs not contain 
X- In both cases. Equations (6.28) and (6.29) indicate that the s = term of the first 
piece and the n = term of the second piece of should be kept. For the piece 
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that is proportional to A^p~^ the exponents are a = 3,p = — 3, x = for the term 

that does not contain x ^'^d a = 3,p = —3, X = 1 for the term that does contain x- 
In this case, Equations (6.28) and (6.29) give that none of the s-terms and only the 
n = term of the second piece of Ef^'^ need to be kept. That gives 



r(o) = e° J]A^Yi^(e,$) 



where, from Equation (6.23) 
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(6.44) 



Substituing the exphcit expressions of the coefficients Ef'^ into Equation (6.44) for 
Bijn results in 



Im 



To - 3M 
2rl{r-o - 2M) 



k,v,X 



X 



N 



(-l)'+"22"+i r(A^ + l) 



n=0 



AIA 



Tq 3Af j^A,m L^-^'™ 

(2n + l)!! r(A^-n + l) V2ro(ro-2M) 3/2 1/2 



A^(ro-3M)5/ ro-3M 



ro (ro 



-iX.m 



-2M) ^'^ro{To-2M) 



+ 



A= 



(r„-3M)2 



\^\^f2r}{ro-2M) 
It is easy to recognize that 
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so the last two terms in the last expression for Bi^ can be combined to give a signifi- 
cantly simpler expression, namely 
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(6.47) 



As is the case for Aim, only the coefficients for which m is even contribute to the sum 
in Equation (6.43). 

The regularization parameter Bj. calculated in [23] can be easily derived from the 
general expression given in Equation (6.47). First, the m = coefficient is considered, 
which also makes A = and u = 0. The first derivative of the factor A|A| with 
respect to r in Equation (6.47) gives zero in the coincidence limit r To- So the only 
term that contributes is the one that is proportional to A. If the spherical harmonic 

y;o(0, $) is again replaced with sj^^ the result is 
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The two finite sums over k and over n can be easily calculated for any vahic of / by 
using Maple. For any value of / that was tried, the result was equal to ^Y^ , but, 



21+1 



again, there is no general proof for that. Substituting this result for the double sum 
into Equation (6.48) gives that 
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which is the same as the result given in [23]. 
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The contribution of T(o) to the derivative of the self-force can now be calculated. 
The second derivative of Equation (6.47) has to be taken, and the m = component 
to be considered. Thus, the term proportional to A vanishes but the second derivative 
of the term A| A| gives 

^(A|A|) = 2sgn(A) (6.50) 
which does give a contribution in the coincidence limit. The result is 
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Using Equation 15.2.13 of [37] and the previously derived result for the hypergeometric 
function Fi (Equation (6.38)), a simple expression for the hypergeometric function F2 
can also be obtained and it is 
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The final expression for the contribution of Bi^n to the second derivative of the scalar 
field is 
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The term of order in Equation (6.22) is considered next. The part of it that 
needs to be kept is 
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For the part that is proportional to p the exponents are a = 0,p = 1, g = 1 for 
the term containing (1/x) and a = 0,p = l.q = 2 for the term containing (1/%^). 
For both these terms, Equations (6.28) and (6.29) indicate that only the s = and 
s = 1 terms of the first piece and none of the values of n of the second piece of 
E^'^ should be kept. For the part that is proportional to A^p~^ the exponents are 
a — 2,p — —1, 5 = for the term not containing x ^-nd a — 2,p — —1, ? = 1 and 
a = 2,p = — 1, g = 2 for the other two terms. In this case, only the s = term should 
be kept in Ef^'^. Consequently 



(6.55) 
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where only the even m's contribute to the above sum and, from Equation (6.23), the 
coefficients Cim are: 
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(6.56) 



After the explicit expressions for the coefficients Ef'^ and for 7^ are substituted into 
the expression for Cim, the result is 



Im 



3M/ 



8ro V r„ 



[ 2V2 / 
l2n + 3r V2 



2M 



EE 



N 



_iy+n22n+| r(A^ + l) 

(2n+l)!! r(Ar-n + l) 



X 



^--^ ^3/2 



+ 



V2ro{ro - 2MY 



2n + 1 fr,, 



2n + 3 



.3M)/z,A,m 
-> V 3, 



3/2 



(6.57) 
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Obtaining the result of [23] for the regularization parameter Cr for the self- 
force is trivial in this case. That is because taking the first derivative of Cim with 
respect to r immediately makes the term with no A dependence vanish and the terms 
proportional to A^ end up having a factor of A which also vanishes in the coincidence 



70 



limit. That means that the 0(e) term does not contibute to the self-force on the 
scalar charge. 

Things arc different for the first derivative of the self-force. Taking the second 
derivative of Equation (6.57) makes the term with no A dependence vanish, but the 
term proportional to does give a contribution which must be calculated. Since the 
coincidence limit is taken, m — and 
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which, after a short calculation, gives 
92 
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The last term of Equation (6.22) that is examined is the term of order e^. The 
part of it that must be considered is 
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In that term, the exponents for A and p are a — l,p — 1, while the exponent of x 
takes the values g = 0, 1, 2, 3 for each term. Equations (6.28) and (6.29) indicate that 
only the s = term of the first piece containing the r-dependence of each Ej'^ needs 
to be kept and the second piece can be disregarded. Then, the O(e^) term can be 
written as 

T(2) = e2^Amllm(©,^) (6.61) 
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where, according to Equation (6.23) 
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Substituting of the coefficients Ef'^ into Equation (6.62) gives 
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Reproducing the result of [23] for Dr can be done by following the same method 
that was used for the regularization parameters Ar and B^.. First, the m = coefficient 
is considered, which makes A = and z/ = 0. The first derivative with respect to r is 
then taken and that simply makes the factor A vanish. Finally the expression of YJo 
for = is used. The result is 
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which is the same as the result of 
used to evaluate the double sum: 



(6.64) 

To derive this final expression. Maple was 
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for different values of Even though a general proof that this equality holds for every 
/ cannot be provided, this equation was verified to hold true for any value of I that 
was tried. 

As far as the contribution of this term to the derivative of the self-force is 
concerned, it can be easily seen to equal zero since there are no terms that can 
survive after the second derivative of Di^ with respect to r is taken. Consequently 



i2 



9- ^ 



= 0. (6.66) 



One final comment should be made about Equations (6.41), (6.53) and (6.59). 
Because these equations involve complicated sums of k and n it may be thought that 
they are not particularly handy for the calculation of the contributions of the cor- 
responding terms to the first derivative of the self-force. Nonetheless, all the sums 
involved are finite and can always be calculated using mathematical software or nu- 
merical code, which renders these equations relatively easy to use. 



CHAPTER 7 
CALCULATION OF THE RETARDED FIELD 



As was explained in Cliapter 4, when calculating the self-force on a scalar charge 
q one needs to know the retarded and the singular part of the scalar field. In order to 
be able to produce results for the self-force, I present here the calculation of the re- 
tarded field generated by a scalar charge moving on a circular orbit in a Schwarzschild 
background, as well as its contribution to the self-force. The method used is similar 
to that presented by Burko in [31] and the calculation was outlined in appendix E 
of [23]. 

7.1 Analytical Work 

The scalar field ■0'^^* at a point p generated by a scalar charge q moving on a 
circular orbit of radius Tq on the equatorial plane of a Schwarzschild background of 
mass M obeys Poisson's equation 

V'V''*(p(a;")) = -47r^ (7.1) 

where x"" are the Schwarzschild coordinates of point p and the source term is given by 
the integral 

-1 J (-9)-'^ ^"^'' 'J'''^ - \m - ^QS{ts - ts{T))dT (7.2) 

^q'±^sie-\)si,-nt.){^y\ 

In this expression p' is the location of the charge q, r is the proper time along the 
circular orbit of the charge q and the relationship between the proper time r and the 
Schwarzschild time ts gives 

dts 1 /„ „s 

ar 3M 

ro 
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Then, by writing Poisson's equation in Scwarzschild coordinates one gets 
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As was mentioned in Chapter 6, both the source and the field can be decomposed 
into spherical harmonics. Specifically, the spherical harmonic decomposition of the 
source can be done using the spherical harmonic decomposition of the 3-dimensional 
S function which is given in [38] 
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where 
and 
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The retarded field can be written as 
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(7.8) 
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where the time dependence is the one implied by the spherical harmonic decomposition 
of the 5 function. 

Substituting the expressions for qim and i/)'^'^^ into Equation (7.4) results in a 
differential equation for the radial part ipimii^) only: 
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For points with r ^ To the right hand side of Equation (7.9) is equal to zero and the 
differential equation can be solved numerically. 
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A variable that is used in the following is r*, which is defined as 
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The behavior of ipim at the two boundaries needs to be considered first. For 
large distances away from the orbit, i.e. in the limit r — > oo, the ipimii^) is expected 
to behave as an outgoing wave, so it is assumed that 
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This expression must satisfy Equation (7.9). The coefficients can be found by 
asserting exactly that. First, the first and second derivatives of ipj^{r) with respect 
to r are calculated. These derivatives are 
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Substituting Equation (7.12) and its two derivatives into Equation (7.9) results to an 
equation for the coefficients a^j. That equation needs to be multiplied by (r — 2MY 
in order to get rid of the (r — 2M) and (r — 2MY factors that appear in some 
denominators. Any factors of (r — 2M) that show up in the numerators can be 
absorbed into the sums of r" so that the final expression contains sums with powers 
of r only. By setting the factor of each r" term equal to zero, a recursion relation for 
the aJs is obtained and it is 



/(/ + 1) - n(n - 1) + UujM 



+ i2M' 



2iujn 
2,n-2)2 

ujn 



a„_i - iM 



(n-l)(2n-3)-/(/ + l) 



O'n-2 



(7.15) 



On-3- 



For this recursion relation it is assumed that Oq = 1 and that a„ = for n < 0. 
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The behavior of il^imif) at the event horizon r = 2M is now considered. There, 
i^imif) should behave as an ingoing wave, so it can be assumed to be 
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(7.16) 



This expression must also be a solution of Equation (7.9). As in the previous case, 
the derivatives are calculated first 
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and 
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Substituting Equation (7.16) and these derivatives into Equation (7.9) gives an equa- 
tion for the coefficients 6„. This equation must first be multiplied by so that the 
factors r^, and r that appear in some denominators vanish. After that, any factors 
of and r that show up in the numerators can be dealt with by the substitutions 



(r - 2Mf + 4M(r - 2M) + 4M^ 
(r - 2M) + 2M 
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so that the final equation contains only sums of powers of (r — 2M). By setting the 
coefficients of each power of (r — 2M) equal to zero, the recursion relation for the 
coefficients b„ is obtained 
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where bo = 1 and 6^ = for n < 0. 

The expression (7.12) with the coefficients (7.15) is used as the starting solution 
for the numerical integration of Equation (7.9) from infinity to the radius of the orbit 
To- This expression is an asymptotic expansion so it is important that the integration 
starts at a value of r which is large enough for the terms of the sum to reach machine 
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accuracy before the series begins to diverge. The result of this first integration is a 
homogeneous sohition of Equation (7.9) for the region r > denoted here by ip^^i''^)- 
Then, the expression (7.16) with the coefficients (7.21) is used as the starting solution 
for the numerical integration of Equation (7.9) from the event horizon r = 2M to the 
radius of the orbit Tq. In this case, the integration starts at a radius r that is close 
enough to 2M for the terms of the sum to reach machine accuracy before the series 
starts to diverge. The result of this integration is a homogeneous solution of Equation 
(7.9) for the region r < To, denoted here as ip^{r). All these solutions are shown in 
Figure (7.1). 



r = 2M 




Figure 7.1: Solutions obtained with the numerical code 

Once the two solutions are found, they have to be matched at the location of the 
particle r = Tq. If ipim{r) stands for the final normalized solution, which results after 
multiplying ipj^{r) and il^i^{r) by the normalization factors Aim and Bim respectively, 
then 

t Bim^^{r), for r>ro 

the factors A^n and Bim can be calculated. The " is used to avoid confusion with the 
regularization parameters Ai^, and Bim calculated in Chapter 6. 
The solution itself should be continuous at r = To- 

Aimi^t,{ro)^Bimi^Z{ro) (7.23) 



but its first derivative with respect to r , also provided by the numerical integration, 
should be discontinuous. The discontinuity is determined by the differential Equation 
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(7.9) when the right-hand side term that contains the 5 function is not set equal to 
zero. Integrating Equation (7.9) from r~ to rj" gives 
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or, taking into account the fact that for all the terms of the left-hand side other 
than the first one the integrand is continuous at r = Tq and the integrals give no 
contribution 
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Equations (7.23) and (7.24) can be easily solved for the factors Aim and and the 
result is 



-^Im — 
Blm — 



Qlr. 



ro - 2M ^H^(r„) '(r,) - V^,- (r„) C(r.) 



ro - 2M ^f^{ro) {ro) - i^Ziro) ^fM 



(7.25) 
(7.26) 



where the primes denote derivatives with respect to r. 

After the retarded field is calculated, its contribution to the self-force along 
the radial direction can be obtained. In fact, every Zm-component as well as every 
Z-component can be calculated by 
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Finally, the total contribution from the retarded field is 
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(7.29) 



1=0 



A note on the angles 6 and (j) is necessary at this point. The usual Schwarzschild 
angles 9 and 4> were used for the calculation of the /-modes of the retarded field in this 



79 



chapter, while the rotated angles © and $ were used for the calculation of the singular 
field and its regularization parameters in Chapter 6. Questions may arise from the 
fact that the two results must be subtracted for the calculation of the self- force. But, 
as was mentioned in Chapter 6, the rotation (6*, 4>) — > (0, $) preserves the spherical 
harmonic index Z so, as long as the summation over m precedes the subtraction of the 
Z-modes, there is no problem with the method. 

7.2 Numerical Code 

The numerical code that was written for the calculation of the radial part of 
the retarded field and its contribution to the self-force is given in Appendix C. The 
language used is C. For the integration of Equation (7.9) the adaptive stepsize Runge- 
Kutta method is used, which is described in [39]. The functions of [39] that were used 
for the integration are: odeintO, rkqsO, rkck(). Also, the function plgndrO 
for the calculation of the Associated Legendre polynomials, also available in [39], was 
used. The details for all these functions are given in [39] and for that reason they are 
not presented here. 

An effort was made to keep the symbols used in the code for the different 
variables identical to the ones used in the analysis of this chapter. In the cases where 
that was not possible, the comments on the code should make the notation clear 
enough for the reader to follow. 



CHAPTER 8 
APPLICATIONS AND CONCLUSIONS 



Some applications of tlie self-force caclulation are presented in this chapter and 
some important conclusions are drawn. The effects of the scalar self-force are pre- 
sented in detail. The case of the gravitational self-force is also discussed. 

8.1 Equations of Motion 

In Chapter 4 it was explained how the self-force on different particles moving 
in curved spacetime can be calculated. Since the ultimate goal is the determination 
of the motion of those particles, it is important to know how the self-force affects the 
equations of motion as well. 

Assume that a point particle that carries a finite charge q (which can be a scalar 
charge, an electric charge or the mass of the particle) is moving in a known background 
spacetime which is characterized by the metric gab- For simplicity also assume that 
no other (external to the particle) scalar, electromagnetic or gravitational fields exist. 
If the effects of the charge q are ignored, the point particle moves on a background 
geodesic. If the geodesic is described by -2'*(s), where s is the proper time and u"' is 
the tangent to the geodesic, the geodesic equation is 

u'^Vau'' = (8.1) 

where Va is the covariant derivative with respect to the background metric gab- The 
geodesic is assumed to be parametrized so that 

u'^Ua = -1. (8.2) 

If the interaction of the particle with its own fields is to be taken into account, 
the self-force T^^^f should be included into the equations of motion. Specifically, if m 
is the mass of the particle, the equation that holds is 

mu^Vau'' = J^.^'eif (8.3) 

where J^^^^^ is given by Equations (4.16), (4.25) and (4.34) for the scalar, electromag- 
netic and gravitational field respectively and is of O(g^). 

8.2 Effects of the Scalar Self-Force 

In general. Equation (8.3) gives the effect of the self-force on each component 
of the 4- velocity of a scalar particle moving in spacetime, when the correct expression 
for the scalar self-force is substituted. If the spherically symmetric Schwarzschild 
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background is considered, it is easily deduced that only 3 components of the 4- velocity 
are affected by the self-force. Indeed, as is explained in [40] , the Schwarzschild metric 
has parity reflection symmetry ^ vr — so if the initial position x"- and the initial 
tangent u°' both lie on the equatorial plane {6 = ^) then the entire path must lie 
on the equatorial plane as well. In addition, due to the spherical symmetry of the 
Schwarzschild background, every path can be brought by a rotation to the equatorial 
plane. Consequently one can only consider equatorial orbits, without imposing any 
real restriction on the path of the particle. The analysis that follows will be presented 
in [41] in more detail. 

Assume that a trajectory of the Schwarzschild background is described by the 
coordinates 

: 0, = T{s),r = R{s), 0^\A- (8-4) 

If E is the energy per unit mass and J the angular momentum per unit mass of the 
particle, it can easily be shown that 

„ dx" r ER ^ J 



ds ^ \(R-2My ' ' R^) 

RR ^^'^^ 

("^' (R-2My °' 



RR 

Ua = gabU 



{R - 2M) 

The normalization condition given in Equation (8.2) for the 4- velocity gives a rela- 
tionship between E, J and i?, specifically 

If the interaction of the scalar particle with its own field is to be taken into account, 
the geodesic equation is modified to 

u^VaUb = qy^il^^. (8.7) 

Calculating the non- vanishing components of Equation (8.7) gives three separate equa- 
tions, namely 

dk M i? - 3M . - 2M ^ , R 

In order to predict the effect of the self-force on circular orbits, one can assume 
that -R = and i? = at some instant. Then Equation (8.10) gives the simple 
expression 

^2 R"^ fM R-2M 



R-3M\R^ R 



-qdr^l^'^) (8.11) 
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which gives the effect of the self-force on the angular momentum per unit mass of the 
scalar particle. 

Also, Equation (8.6) gives 

, / 2M\ / J2\ ^ ^ 



R J\ E?)' 
Combining Equations (8.11) and (8.12) results in 

which shows the effect of the self- force on the energy of the scalar particle. All the 
results mentioned up to this point are exact. 

Since the effect of the scalar charge q can be assumed to be small, the approxi- 
mation that (Rdri/j^) « 1 can be made, and Equation (8.13) for the energy E can 
be written as 

E = —^Z2^={1 - IqRdri^^). (8.14) 
The orbital frequency is defined as 

from Equation (8.5). Substituting Equations (8.11) and (8.14) into Equation (8.15), 
and using again the approximation (Rdrip^) << 1, a simple expression for the effect 
of the self-force on the orbital frequency is obtained: 



0' = 




2 M 



As expected, if the interaction of the particle with its own scalar field is ignored, the 
well-known result Q — y/MR~^ is recovered. 



8.2.1 Calculation of the Self- Force 

A sample self-force calculation is presented in this section. The scalar charge is 
assumed to be g = 1 and to be moving on a circular orbit of radius R — 10 M in a 
Schwarzschild background of M = 1. 

The Zm-components ipimi'r) of the retarded field generated by the scalar charge 
and their contributions Tlf to the Z-multipole modes of the retarded part of the self- 
force arc first calculated for / = 0, 1, . . . ,40, using the numerical code described in 
Chapter 7 and given in Appendix C. Then, the non-zero regularization parameters A^., 
Br and Dj. are used ^ to subtract off the contributions of the singular-field self-force 
up to 0(/^^), according to Equation (6.12). These results are shown in Figure (8.1). 

^Remember that it was proven that C^. = for circular Schwarzschild orbits, in Chapter 6. 
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In that figure, the curve marked Fi is the contribution J^f' as a function of I. The 
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Figure 8.1: Scalar Self- Force 

curves marked A, B and D show the terms {J^lf — J-'fr) functions of where J-'f^. 
successively includes the contributions from the regularization parameters Ar, Br and 
Dr respectively. It is clear that the series for the self-force converges when the terms 
containing the regularization parameter B^ are included. However, the convergence 
becomes much faster the terms containing the parameter Dj. are also included. 

As is rigorously explained in [23], the contributions from the terms in Equa- 
tion (6.12) are successively of order Specifically, the 0{l^^) term is 

associated with the parameter E^, the 0{l~^) term is associated with the parameter 
E^ and, in general, the 0{l~'^^) term is associated with the parameter E"^"^. The 
parameters E^ are /-independent [23] and can be determined as follows. The contri- 
bution from each term is proven [23] to be 
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(8.18) 



Assume that the residual {TJ^^ ~ Tf^) as a function of resulting after the subtraction 

is calculated. The coefficients E^ can be 



of the contribution of Ar 



Br and D^, 
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determined numerically by fitting that residual, to a linear combination of terms of 
the form (8.17) and successively removing their contributions. 

The curves marked E^,. . . in Figure (8.1) are the residuals resulting after 
the numerical fit of the parameters E^, . . . respectively. It was noticed that fitting 
more than four of the coefficients E^ did not improve the residual, which had already 
reached machine accuracy by the fourth fitting. 



8.2.2 Change in Orbital Frequency 

It was shown in Section (8.2) that the orbital frequency of the circular orbit is 
affected by the self-force. Using Equation (8.16) the difference in frequency can be 
calculated; it is 

1 R{R - 3M) 
2 



U9) 



n n 2 M 

The effect of the self-force on the orbital frequency is shown in Figure (8.2), where 
the ratio (AQ/Q) is plotted as a function of (R/M). The orbital frequency Q before 
the self-force effects are included is also plotted for comparison. 




R/M 



Figure 8.2: Change in Orbital Frequency 



8.3 Effects of the Gravitational Self-Force 

Assume that a particle of mass m is moving in a background spacetime described 
by the metric gab, with no other extremal fields present. Equation (8.3) gives the 
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worldline of the particle that includes the self-force effect if ^g"if is substituted from 
Equation (4.34). Clearly, given the initial position and 4- velocity of the particle, 
that worldline is determined exclusively by the background metric Qah and the metric 
perturbation h^^^ induced by the particle. Both of these contributions are solutions of 
the homogeneous Einstein equations. 

An observer making local measurements of the metric in the vicinity of the 
particle measures two contributions to it. The first contribution comes from the 
background metric Qat combined with the /i^^^-part of the actual metric perturbation 
generated by the particle. Since both these terms obey the homogeneous Einstein 
equations, the observer sees no source for either one, so no local measurement can allow 
him to distinguish /i^^ from the background metric gab- The second contribution is /i^^, 
a metric perturbation also generated by the particle which, according to the observer, 
is sourced by the particle. The local observer sees the particle move in a background 
geodesic of the metric {gab + h^b)- Because he only makes local measurements, he 
observes no radiation generated by the moving particle and consequently no effect 
that he could describe as radiation reaction [22] . 



APPENDIX A 
GRTENSOR CODE FOR THE SINGULAR FIELDS 



restart ; 

readlib(grii) ; 
grtensor () ; 

qload(thz) ; 

e := epsilon; 

ord := [t, X, y, z, Bxx, Byy] ; 

EpsEqn := {e"10=0, e'^9=0, e'"8=0, e''7=0, e"6=0, e''5=0, e'^4=0, e"3=0, 

e~2=0}; 

Traces := {Ezz=-Exx-Eyy , Bzz=-Bxx-Byy} ; 



# The Ricci Tensor Should Be Zero # 

grcalc (detg) ; 

grmapCdetg, series, 'x', e, 4); 
grmapCdetg, convert, 'x' , polynom) ; 

grmapCdetg, expand, 'x'); 
grmapCdetg, subs, EpsEqn, 'x'); 

grmapCdetg, collect, 'x' , ord, distributed, factor); 
grmapCdetg, subs. Traces, 'x'); 

grmapCdetg, collect, 'x' , ord, distributed, f actor) ; 
grdisplayCdetg) ; 

grcalcCgCup,up)) ; 

grmapCgCup,up) , series, 'x', e, 4); 
grmapCgCup,up) , convert, 'x' , polynom); 
grmap Cg Cup , up) , expand , ' x ; 
grmapCgCup,up) , subs, EpsEqn, 'x'); 

grmap Cg Cup, up) , collect, 'x' , ord, distributed, factor); 
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grmap(g(up,up) , subs. Traces, 'x'); 

grmap(g(up,up) , collect, 'x' , ord, distributed, factor); 
grdisplay(g(up,up)) ; 

grcalc (Chr (dn , dn , dn) ) ; 

grmapCChr (dn,dn,dn) , series, 'x', e, 4); 
grmap (Chr (dn , dn , dn) , convert , ' x ' , polynom) ; 
grmap(Chr(dn,dn,dn) , expand, 'x'); 
grmap (Chr (dn,dn,dn) , subs, EpsEqn, 'x'); 

grmap (Chr (dn,dn,dn) , collect, 'x', ord, distributed, factor 
grmap (Chr (dn,dn,dn) , subs. Traces, 'x'); 

grmap (Chr (dn,dn,dn) , collect, 'x', ord, distributed, factor 
grdisplay (Chr (dn , dn , dn) ) ; 

grcalc (Chr (dn,dn, up) ) ; 

grmap (Chr (dn,dn, up) , series, 'x', e, 4); 

grmap (Chr (dn,dn, up) , convert, 'x' , polynom); 
grmap (Chr (dn,dn, up) , expand, 'x'); 

grmap (Chr (dn,dn, up) , subs, EpsEqn, 'x'); 

grmap (Chr (dn,dn, up) , collect, 'x', ord, distributed, factor 

grmap (Chr (dn,dn, up) , subs. Traces, 'x'); 

grmap (Chr (dn,dn, up) , collect, 'x' , ord, distributed, factor 
grdisplay (Chr (dn,dn, up) ) ; 

grcalc (Chr (dn , up , up) ) ; 

grmap (Chr (dn, up, up) , series, 'x', e, 4); 

grmap (Chr (dn, up, up) , convert, ^x' , polynom); 
grmap (Chr (dn , up , up) , expand , 'x') ; 

grmap (Chr (dn, up, up) , subs, EpsEqn, 'x'); 

grmap (Chr (dn, up, up) , collect, 'x' , ord, distributed, factor 

grmap (Chr (dn , up , up) , subs , Traces , ' x ' ) ; 

grmap (Chr (dn, up, up) , collect, 'x', ord, distributed, factor 
grdisplay (Chr (dn , up , up) ) ; 

grcalc (Chr (dn , up , dn) ) ; 

grmap (Chr (dn, up, dn) , series, 'x', e, 4); 

grmap ( Chr (dn, up, dn) , convert, 'x', polynom); 
grmap (Chr (dn, up, dn) , expand, 'x'); 

grmap (Chr (dn, up, dn) , subs, EpsEqn, 'x'); 

grmap ( Chr (dn, up, dn) , collect, 'x' , ord, distributed, factor 

grmap (Chr (dn, up, dn) , subs. Traces, 'x'); 

grmap (Chr (dn, up, dn) , collect, 'x' , ord, distributed, factor 
grdisplay (Chr (dn , up , dn) ) ; 
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grcalc(R(dn,dn,dn,dn)) ; 

grmap(R(dn,dn,dn,dn) , series, 'x', e, 4); 
grmap(R(dn,dn,dn,dn) , convert, 'x', polynom) ; 
grmap(R(drL,dn,dn,dn) , expand, 'x'); 
grmap(R(dn,dn,dn,dn) , subs, EpsEqn, 'x'); 

grmap(R(dn,dn,dn,dn) , collect, 'x' , ord, distributed, factor); 

grmap(R(dn,dn,dn,dn) , subs. Traces, 'x'); 

grmap(R(dn,dn,dn,dn) , collect, 'x', ord, distributed, factor); 
grdisplay (R(dn,dn,dn,dn) ) ; 

grcalc(R(up,up,dn,dn) ) ; 

grmap(R(up,up,dn,dn) , series, 'x' , e, 4); 

grmap(R(up,up,dn,dn) , convert, 'x' , polynom); 
grmap(R(up,up,dn,dn) , expand, 'x'); 

grmap(R(up,up,dn,dn) , subs, EpsEqn, 'x'); 

grmap(R(up,up,dn,dn) , collect, 'x' , ord, distributed, factor); 

grmap(R(up,up,dn,dn) , subs. Traces, 'x'); 

grmap(R(up,up,dn,dn) , collect, 'x' , ord, distributed, factor); 
grdisplay (R (up, up, dn,dn) ) ; 

grcalc(R(up,dn,dn,dn) ) ; 

grmap(R(up,dn,dn,dn) , series, 'x' , e, 4); 

grmap(R(up,dn,dn,dn) , convert, 'x' , polynom); 
grmap(R(up,dn,dn,dn) , expand, 'x'); 

grmap(R(up,dn,dn,dn) , subs, EpsEqn, 'x'); 

grmap(R(up,dn,dn,dn) , collect, 'x' , ord, distributed, factor); 

grmap(R(up,dn,dn,dn) , subs. Traces, ^x'); 

grmap(R(up,dn,dn,dn) , collect, 'x', ord, distributed, factor); 
grdisplay (R (up, dn,dn,dn) ) ; 

grcalc(R(dn,up,dn,up) ) ; 

grmap(R(dn,up,dn,up) , series, 'x' , e, 4); 

grmap(R(dn,up,dn,up) , convert, ^x' , polynom); 
grmap(R(dn,up,dn,up) , expand, 'x'); 

grmap(R(dn,up,dn,up) , subs, EpsEqn, 'x'); 

grmap(R(dn,up,dn,up) , collect, 'x' , ord, distributed, factor); 

grmap(R(dn,up,dn,up) , subs. Traces, ^x'); 

grmap(R(dn,up,dn,up) , collect, 'x' , ord, distributed, factor); 
grdisplay (R(dn, up, dn, up) ) ; 

grcalc(R(dn,dn)) ; 

grmap(R(dn, dn) , series, 'x', e, 4); 
grmap(R(dn,dn) , convert, 'x' , polynom); 

grmap(R(dn,dn) , expand, 'x'); grmap(R(dn,dn) , subs, EpsEqn, 'x'); 
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grmap(R(dn,dn) , collect, 'x' , ord, distributed, factor); 

grmap(R(dn,dn) , subs, Traces, 'x'); 

grmap(R(dn,dn) , collect, 'x' , ord, distributed, factor); 
grdisplay (R(dn,dn) ) ; 



# Useful Quantities # 

grdef ('etadn{(a b)} := -kdelta{''$t a}*kdelta{'"$t b} 
+kdelta-["$x a}*kdelta{~$x b} + kdelta{''$y a}*kdelta{'"$y b} 
+kdelta{"$z a}*kdelta{~$z b}'); 
grcalc(etadn(dn,dn) ) ; 
grdisplay (etadn(dn,dn) ) ; 

grdef ('etaup{("a "b)} := -kdelta{$t ~a}*kdelta{$t ~b} 
+kdelta-[$x "a}*kdelta{$x "b} + kdelta{$y ~a}*kdelta{$y "b} 
+kdelta{$z "a}*kdelta{$z ~b}'); 
grcalc(etaup(up,up)) ; 
grdisplay (etaup (up, up) ) ; 

grdef ('eps{[a be]} := LevCS{d a b c}*kdelta{$t "d}'); 
grcalc(eps(dn,dn,dn)) ; 
grdisplay (eps(dn,dn,dn)) ; 

grdef ('EE-[(a b)} := kdelta{"$x a}*kdelta{~$x b}*e*Exx 
+ kdelta{"$y a}*kdelta{~$y b}*e*Eyy 
+ kdelta{''$z a}*kdelta{'^$z b}*e*Ezz 
+ kdelta{''$x a}*kdelta{'^$y b}*e*Exy 

+ kdelta{~$x a}*kdelta{"$z b}*e*Exz 
+ kdelta-["$y a}*kdelta{~$z b}*e*Eyz'); 
grcalc (EE (dn , dn) ) ; 
grdisplay (EE (dn , dn) ) ; 

grdef ('BB-[(a b)} := kdelta{"$x a}*kdelta{~$x b}*e*Bxx 

+ kdelta-["$y a}*kdelta{~$y b}*e*Byy 

+ kdelta{"$z a}*kdelta{~$z b}*e*Bzz 

+ kdelta{~$x a}*kdelta{~$y b}*e*Bxy 

+ kdelta{''$x a}*kdelta{~$z b}*e*Bxz 

+ kdelta{"$y a}*kdelta{'^$z b}*e*Byz'); 

grcalc (BB(dn,dn)) ; 

grdisplay (BB (dn , dn) ) ; 
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grdef ('hO-[(b c)} : = 

kdelta{~$t b}*kdelta{"$x c}*2*A*mu*y/((x~2+y"2+z"2)"(3/2)) 
+ kdelta{"$t b}*kdelta{''$y c}*(-2)*A*inu*x/((x~2+y~2+z''2)''(3/2)) ') ; 
grcalc(hO(dn,dn) ) ; 
grdisplay (hO (dn , dn) ) ; 

grdef ( ' epsx{ [a b]} := eps{c a b}*kdelta{$x "c}'); 
grcalc (epsx (dn , dn) ) ; 
grdisplay (epsx (dn, dn) ) ; 

grdef ( 'epsy{ [a b]> := eps{c a b}*kdelta{$y "c}'); 
grcalc (epsy (dn , dn) ) ; 
grdisplay (epsy (dn, dn) ) ; 

grdef ( 'epsz-[ [a b] } := eps{c a b}*kdelta{$z "c}'); 
grcalc (epsz(dn,dn)) ; 
grdisplay (epsz (dn, dn) ) ; 



# Chech What The Equations Look Like: # 

grdef ('hltest{(a b)} := kdelta{"$t a}*kdelta{'"$t b}*e*hltt 

+ kdelta-["$t a}*kdelta{~$x b}*e*hltx 

+ kdelta-["$t a}*kdelta{~$y b}*e*hlty 

+ kdelta{"$t a}*kdelta{~$z b}*e*hltz 

+ kdelta{''$x a}*kdelta{'^$x b}*e*hlxx 

+ kdelta{''$x a}*kdelta{''$y b}*e*hlxy 

+ kdelta-["$x a}*kdelta{~$z b}*e*hlxz 

+ kdelta-["$y a}*kdelta{~$y b}*e*hlyy 

+ kdelta{"$y a}*kdelta{~$z b}*e*hlyz 

+ kdelta{''$z a}*kdelta{'^$z b}*e*hlzz'); 

grcalc (hltest (dn , dn) ) ; 

grdisplay (hltest(dn,dn)) ; 

grdef ('htest{(a b)} := hO{a b} + hltest{a b}' ); 
grcalc (htest (dn , dn) ) ; 
grdisplay (htest (dn , dn) ) ; 

grdef ( 'htesttrace := htest{"c c}'); 
grcalc (htesttrace) ; 
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grmapChtesttrace, series, 'x' , e, 4); 

grmapChtesttrace , convert, 'x', polynom) ; 

grmapChtesttrace, expand, 'x'); grmapChtesttrace, subs, EpsEqn, 'x'); 
grmapChtesttrace, collect, 'x^ , ord, distributed, factor); 
grmapChtesttrace, subs. Traces, 'x'); 

grmapChtesttrace, collect, 'x' , ord, distributed, factor); 
grdisplayChtesttrace) ; 

grdef C'hbartest{Ca b)} := htest{a b} - Cl/2)*g{a b}*htesttrace ' ) ; 
grcalc Chbart est Cdn , dn) ) ; 

grmapChbartestCdn,dn) , series, 'x' , e, 4); 
grmapChbartestCdn,dn) , convert, 'x', polynom); 
grmap Chbartest Cdn , dn) , expand , 'x') ; 
grmapChbartest Cdn,dn) , subs, EpsEqn, 'x'); 

grmap Chbartest Cdn, dn) , collect, 'x' , ord, distributed, factor); 
grmap Chbartest Cdn, dn) , subs. Traces, 'x'); 

grmap Chbartest Cdn, dn) , collect, 'x', ord, distributed, factor); 
grdisplay Chbartest Cdn , dn) ) ; 



# Left-Hand Side Of The Equation # 

grcalc Chbartest Cdn , dn , cdn) ) ; 

grmap Chbartest Cdn , dn , cdn) , series , ' x ' , e , 4) ; 

grmap Chbartest Cdn , dn , cdn) , convert , ' x ' , polynom) ; 

grmap Chbartest Cdn, dn, cdn) , expand, 'x'); 

grmap Chbartest Cdn, dn, cdn) , subs, EpsEqn, 'x'); 

grmapChbartest Cdn, dn, cdn) , collect, 'x', ord, distributed, factor); 

grmap Chbartest Cdn , dn , cdn) , subs , Traces , ' x ' ) ; 

grmap Chbartest Cdn, dn, cdn) , collect, 'x', ord, distributed, factor); 
grdisplay Chbartest Cdn , dn , cdn) ) ; 

grdef C 'Dhbartest{a b "c} := hbartest{a b ;"c}'); 
grcalc CDhbartest Cdn , dn , up) ) ; 

grmap CDhbart est Cdn, dn, up) , series, 'x', e, 4); 
grmap CDhbartest Cdn, dn, up) , convert, 'x', polynom); 
grmap CDhbartest Cdn, dn, up) , expand, 'x'); 
grmap CDhbartest Cdn, dn, up) , subs, EpsEqn, 'x'); 

grmap CDhbartest Cdn, dn, up) , collect, 'x' , ord, distributed, factor); 
grmap CDhbartest Cdn , dn , up) , subs , Traces , ' x ' ) ; 

grmap CDhbartest Cdn, dn, up) , collect, 'x', ord, distributed, factor); 
grdisplay CDhbartest Cdn , dn , up) ) ; 
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grdef ('DDhbartest{a b} := Dhbartest{a b "c ;c}'); 
grcalc(DDhbartest(dn,dn)) ; 

grmap(DDhbartest(dn,dn) , series, 'x', e, 4); 
grmap(DDhbartest (dn,dn) , convert, 'x' , polynom) ; 
grmap (DDhbartest (dn , dn) , expand , ' x ' ) ; 
grmap(DDhbartest(dn,dn) , subs, EpsEqn, 'x'); 

grmap (DDhbartest (dn, dn) , collect, 'x', ord, distributed, factor); 
grmap (DDhbarte St ( dn, dn) , subs. Traces, 'x'); 

grmap (DDhbartest (dn,dn) , collect, 'x', ord, distributed, factor); 
grdisplay (DDhbartest (dn , dn) ) ; 

grdef ('Rhbartest{a b} := 2*R{a "c b "d}*hbartest{c d}'); 
grcalc (Rhbartest (dn , dn) ) ; 

grmap (Rhbartest (dn,dn) , series, 'x', e, 4); 
grmap (Rhbartest (dn,dn) , convert, 'x' , polynom); 
grmap (Rhbcirt est (dn , dn) , expand , ' x ' ) ; 
grmap (Rhbartest (dn,dn) , subs, EpsEqn, 'x'); 

grmap (Rhbartest (dn, dn) , collect, 'x', ord, distributed, factor); 
grmap (Rhbartest (dn,dn) , subs. Traces, 'x'); 

grmap (Rhbartest (dn,dn) , collect, 'x' , ord, distributed, factor); 
grdisplay (Rhbartest (dn , dn) ) ; 

grdef ('LHStest{(a b)} := DDhbartest{a b} + Rhbartest{a b}'); 
grcalc (LHStest(dn,dn)) ; 

grmap (LHStest (dn,dn) , series, 'x' , e, 4); 
grmap (LHStest (dn,dn) , convert, 'x' , polynom); 
grmap (LHStest (dn,dn) , expand, 'x'); 
grmap (LHStest (dn, dn) , subs, EpsEqn, 'x'); 

grmap (LHStest (dn,dn) , collect, 'x' , ord, distributed, factor); 
grmap (LHStest (dn,dn) , subs. Traces, 'x'); 

grmap (LHStest (dn,dn) , collect, 'x', ord, distributed, factor); 
grdisplay (LHStest (dn , dn) ) ; 



# Assume A Specific Form For The Solution # 

grdef ('AngM{~a} := [0, 0, 0, A]'); 
grcalc (AngM (up)) ; grdisplay (AngM (up) ) ; 



grdef ('hl-[(a b)} : = 

kdelta{~$t a}*kdelta{"$t b}* (mu/ ( (x~2+y"2+z"2) ~ (3/2) ) ) * 
(att*etadn{i l}*AngM{"l}*x{"i}*BB{j k}*x{" j}*x{"k} 
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+ btt*AngM{~i}*BB{i j}*x{~j}*(x"2+y''2+z''2)) 
+kdelta-["$t a}*kdelta{~$x b}* (mu/( (x"2+y"2+z"2) ~ (3/2) ) )* 

(cx*epsx{i j}*AngM{~i}*etaup{~j "n}*EE{n l}*x{"l}* (x~2+y~2+z~2) 
+ dx*epsx{i j}*AngM{"i}*x{~j}*EE{l k}*x{~l}*x{~k} 
+ fx*epsx{i j}*AngM{'"l}*etaup{"i ''n}*EE{n l}*x{''j}*(x''2+y"2+z"2) 
+ gx*epsx{i j}*etadn{in l}*AngM{"m}*etaup{"i ''n}*EE{n k} 

*x{~j}*x{~l}*x{~k} 
+ px*eps-[i j l}*AngM{"i}*etaup-["j ~n}*EE-[$x n>*x{~l}*(x"2+y"2+z"2) 
+ qx*eps{i j l}*AngM{"i}*etaup{" j ~n}*EE{n k> 
*x{'^d}*kdelta{''$x d}*x{''l}*x{'^k}) 
+kdelta{''$t a}*kdelta{"$y b}* (mu/((x'"2+y''2+z''2)''(3/2)))* 

(cy*epsy{i j}*AngM{"i}*etaup{~ j ~n}*EE{n l}*x{'"l}*(x"2+y''2+z''2) 
+ dy*epsy{i j}*AngM{"i}*x{~ j}*EE-[l k}*x{~l}*x{~k} 
+ fy*epsy{i j}*AngM{"l}*etaup{~i "n}*EE{n l}*x{" j}* (x~2+y~2+z~2) 
+ gy*epsy{i j}*etadn{in l}*AngM{'"m}*etaup{~i ''n}*EE{n k} 

*x{" j }*x{'"l}*x{''k} 
+ py*eps-[i j l}*AngM{~i}*etaup{~j ~n}*EE-[$y n>*x{'^l}*(x''2+y''2+z''2) 
+ qy*eps-[i j l}*AngM{~i>*etaup-[" j ~n}*EE{n k} 
*x{~d}*kdelta{"$y d}*x{"l}*x{~k}) 
+kdelta{''$t a}*kdelta{'^$z b}* (mu/((x''2+y''2+z''2)'^(3/2)))* 

(cz*epsz{i j}*AngM{''i}*etaup{"j ~n}*EE{n l}*x{'"l}*(x"2+y''2+z''2) 
+ dz*epsz{i j}*AngM{"i}*x{"j>*EE{l k}*x{~l}*x{~k} 
+ fz*epsz{i j}*AngM{~l}*etaup{~i "n}*EE{n l}*x{" j}* (x~2+y~2+z~2) 
+ gz*epsz{i j}*etadn{m l>*AngM{~m}*etaup{~i "n}*EE{n k} 

*x{" j }*x{~l}*x{~k} 
+ pz*eps{i j l}*AngM{~i}*etaup{~j "n}*EE{$z n}*x{~l}*(x''2+y~2+z''2) 
+ qz*eps-[i j l}*AngM{~i}*etaup{" j ~n}*EE{n k} 
*x{~d}*kdelta{"$z d}*x{"l}*x{~k}) 
+kdelta{"$x a}*kdelta{"$x b}* (mu/( (x"2+y"2+z"2) ~ (3/2)) )* 
(cxx*etadn{i $x}*AngM{'"i}*BB{$x k}*x{~k}*(x'^2+y'^2+z''2) 
+ dxx*etadn{i $x}*AngM{'"i}*BB{k l}*x{~k}*etadn{$x j}*x{'" j}*x{"l} 
+ fxx*etadn{k i}*AngM{~i>*BB{$x $x}*x{~k}* (x"2+y"2+z~2) 
+ gxx*AngM-["k}*BB{$x k}*etadn{i $x}*x-["i}* (x"2+y"2+z~2) 
+ pxx*etadn{i k}*AngM{~i}*BB{$x l}*etadn{j $x}*x{" j}*x{"k}*x{~l} 
+ qxx*AngM{''k}*BB{k l}*etadn{i $x}*x{''i}*etadn{j $x}*x{'' j}*x{"l} 
+ vxx*etadn{$x $x}*etadn{i l}*AngM{''l}*BB{j k}*x{~i}*x{"j}*x{''k} 
+ uxx*etadn{$x $x}*AngM{~i}*BB{i j}*x{~ j}* (x~2+y~2+z"2) ) 
+kdelta-["$x a>*kdelta{"$y b}*(mu/((x"2+y~2+z"2) " (3/2)) )* 
(axy*AngM{~i}*kdelta{i "$y}*BB{$x k}*x{~k}* (x~2+y"2+z"2) 
+ bxy*AngM{~k}*kdelta{k ~$y}*BB{i j}*x*x{"i}*x{" j} 
+ cxy*etadn{i $x}*AngM{'"i}*BB{$y k}*x{"k}*(x'^2+y'^2+z''2) 
+ dxy*etadn{i $x}*AngM{"i}*BB-[k l}*x{~k}*etadn{$y j}*x{" j}*x{"l} 
+ fxy*etadn{k i}*AngM-["i}*BB{$x $y}*x-["k}* (x"2+y"2+z~2) 
+ gxy*AngM{"k}*BB{$x k}*etadn{i $y}*x{"i}* (x"2+y"2+z~2) 
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+ pxy*etadn{i k}*AngM{~i}*BB{$x l}*etadn{j $y}*x{'' j}*x{"k}*x{~l} 
+ qxy*AngM-["k}*BB{k l}*etadn-[i $x}*x{~i}*etadn{j $y}*x{" j}*x{"l} 
+ sxy*etadn{i k}*AngM{"k}*BB{$y j}*x*x{"i}*x{~ j} 
+ txy*AngM{"j}*BB{$y j}*x* (x~2+y~2+z~2) ) 

+kdelta{''$x a}*kdelta{~$z b}*(mu/((x''2+y''2+z''2)~(3/2)))* 
(axz*AngM{'^i}*kdelta{i '■$z}*BB{$x k}*x{'"k}*(x'"2+y'"2+z''2) 
+ bxz*AngM-["k}*kdelta{k ~$z}*BB-[i j}*x*x{"i}*x{" j} 
+ cxz*etadn{i $x}*AngM{~i}*BB-[$z k}*x{~k}* (x"2+y"2+z"2) 
+ dxz*etadn{i $x}*AngM{~i}*BB{k l}*x{~k}*etadn{$z j}*x{" j}*x{~l} 
+ fxz*etadn{k i}*AngM{'^i}*BB{$x $z}*x{''k}*(x'"2+y'"2+z'^2) 
+ gxz*AngM{''k}*BB{$x k}*etadn{i $z}*x{~i}*(x'"2+y''2+z'^2) 
+ pxz*etadn{i k}*AngM{~i}*BB{$x l}*etadn{j $z}*x{^ j}*x{^k}*x{"l} 
+ qxz*AngM-C"k}*BB{k l}*etadn-[i $x}*x{~i>*etadn{j $z}*x{^ j}*x{"l} 
+ sxz*etadn{i k}*AngM{"k}*BB{$z j}*x*x{"i}*x{" j} 
+ txz*AngM{~j}*BB{$z j}*x* (x"2+y"2+z'^2) ) 

+kdelta{'^$y a}*kdelta{~$y b}*(mu/((x'^2+y'"2+z''2)"(3/2)))* 
(cyy*etadn-[i $y}*AngM{"i}*BB-[$y k}*x{~k}* (x~2+y~2+z"2) 
+ dyy*etadn{i $y}*AngM{"i}*BB-[k l}*x{~k}*etadn{$y j}*x{" j}*x{~l} 
+ fyy*etadn{k i}*AngM{~i}*BB{$y $y}*x-["k}* (x"2+y"2+z~2) 
+ gyy*AngM-["k}*BB{$y k}*etadn{i $y}*x{''i}*(x'"2+y'"2+z'^2) 
+ pyy*etadn{i k}*AngM{''i}*BB{$y l}*etadn{j $y}*x{'" j}*x{'"k}*x{~l} 
+ qyy*AngM-["k}*BB{k l}*etadn{i $y}*x{~i}*etadn{j $y}*x{" j}*x{~l} 
+ vyy*etadn{$y $y}*etadn{i l}*AngM{"l}*BB{j k}*x{"i>*x-[~ j}*x{"k} 
+ uyy*etadn{$y $y}*AngM{~i}*BB{i j}*x{~j}*(x~2+y~2+z"2) ) 

+kdelta{~$y a}*kdelta{~$z b}*(mu/((x~2+y''2+z''2)~(3/2)))* 
(ayz*AngM{~i}*kdelta{i '■$z}*BB{$y k}*x{''k}*(x'"2+y'"2+z''2) 
+ byz*AngM-["k}*kdelta{k ~$z}*BB-[i j}*y*x{"i}*x{" j} 
+ cyz*etadn{i $y}*AngM{~i}*BB-[$z k}*x{~k}* (x"2+y~2+z"2) 
+ dyz*etadn{i $y}*AngM{"i}*BB{k l}*x{~k}*etadn{$z j}*x{" j}*x{~l} 
+ fyz*etadn{k i}*AngM{^i}*BB{$y $z}*x{~k}*(x'"2+y'"2+z'^2) 
+ gyz*AngM{''k}*BB{$y k}*etadn{i $z}*x{''i}*(x''2+y''2+z~2) 
+ pyz*etadn{i k}*AngM{~i}*BB{$y l}*etadn{j $z}*x{" j}*x{"k}*x{~l} 
+ qyz*AngM-["k}*BB{k l}*etadn-[i $y}*x{~i}*etadn{j $z}*x{" j}*x{~l} 
+ syz*etadn{i k}*AngM{~k}*BB{$z j}*y*x{~i}*x{" j} 
+ tyz*AngM{"j}*BB{$z j}*y*(x"2+y"2+z"2)) 

+kdelta{''$z a}*kdelta{'^$z b}*(mu/((x''2+y''2+z''2)~(3/2)))* 
(czz*etadn-[i $z}*AngM{~i}*BB-[$z k}*x{~k}* (x~2+y~2+z~2) 
+ dzz*etadn{i $z}*AngM{~i}*BB-[k l}*x{~k}*etadn{$z j}*x{" j}*x{~l} 
+ fzz*etadn{k i}*AngM{"i}*BB{$z $z}*x{"k}* (x"2+y"2+z~2) 
+ gzz*AngM-["k}*BB{$z k}*etadn{i $z}*x{"i}* (x"2+y"2+z~2) 
+ pzz*etadn{i k}*AngM{'^i}*BB{$z l}*etadn{j $z}*x{'' j}*x{''k}*x{'^l} 
+ qzz*AngM-["k}*BB{k l>*etadn-[i $z}*x{~i}*etadn{j $z}*x{" j}*x{~l} 
+ vzz*etadn{$z $z}*etadn{i l}*AngM-["l>*BB-[j k}*x{"i}*x-[" j}*x{~k} 
+ uzz*etadn{$z $z}*AngM{~i}*BB{i j}*x{~j}*(x~2+y~2+z"2) ) ' ); 



grcalc (hi (dn , dn) ) ; 



grdef ('hh{(a b)} := hO{a b} + hl{a b}'); 
grcalc (hh (dn , dn) ) ; 

grdef ( 'htrace := hh{~c c}'); 
grcalc (htrace) ; 

grmap (htrace , series, 'x', e, 4); 
grmap (htrace , convert, 'x' , polynom) ; 
grmap (htrace , expand, 'x'); 
grmap (htrace , subs, EpsEqn, 'x'); 

grmap(htrace , collect, 'x', ord, distributed, factor); 
grmap (htrace , subs. Traces, 'x'); 

grmap (htrace , collect, 'x' , ord, distributed, factor); 
grdisplay (htrace) ; 

grdef ('hbar-[(a b)} := hh{a b} - (l/2)*g{a b}*htrace'); 

grcalc (hbar(dn,dn)) ; 

grmap (hbar (dn, dn) , series, 'x' , e, 4); 

grmap (hbar (dn, dn) , convert, 'x' , polynom); 

grmap (hbar (dn , dn) , expand , ' x ' ) ; 

grmap(hbar (dn,dn) , subs, EpsEqn, 'x'); 

grmap(hbar (dn,dn) , collect, 'x' , ord, distributed, factor); 
grmap (hbar (dn, dn) , subs. Traces, 'x'); 

grmap (hbar (dn, dn) , collect, 'x' , ord, distributed, factor); 
grdisplay (hbar (dn , dn) ) ; 

grcalc (hbar (dn , dn , cdn) ) ; 

grmap (hbar (dn, dn, cdn) , series, 'x', e, 4); 

grmap (hbar (dn,dn, cdn) , convert, 'x', polynom); 
grmap (hbar (dn , dn , cdn) , expand , ' x ' ) ; 

grmap(hbar (dn,dn, cdn) , subs, EpsEqn, 'x'); 

grmap (hbar (dn,dn, cdn) , collect, 'x', ord, distributed, factor 

grmap (hbar (dn,dn, cdn) , subs. Traces, 'x'); 

grmap (hbar (dn,dn, cdn) , collect, 'x', ord, distributed, factor 
grdisplay (hbar (dn , dn , cdn) ) ; 

grdef ( 'Dhbar{a b ~c} := hbar{a b ;"c}'); 
grcalc (Dhbar (dn , dn , up) ) ; 

grmap (Dhbar(dn,dn, up) , series, 'x' , e, 4); 
grmap (Dhbar (dn,dn, up) , convert, 'x', polynom); 
grmap (Dhbar (dn , dn , up) , expand , ' x ' ) ; 
grmap (Dhbar (dn,dn, up) , subs, EpsEqn, 'x'); 



grmap(Dhbar (dn,dn,up) , collect, 'x', ord, distributed, factor 

grmap(Dhbar(dn,dn,up) , subs, Traces, 'x'); 

grmap(Dhbar(dn,dn,up) , collect, 'x', ord, distributed, factor 
grdisplay (Dhbar (dn , dn , up) ) ; 

grdef ('DDhbar{a b} := Dhbar{a b ~c ;c}'); 
grcalc(DDhbar(dn,dn)) ; 

grmap(DDhbar (dn,dn) , series, 'x', e, 4); 
grmap(DDhbar (dn,dn) , convert, ''x.' , polynom) ; 
grmap (DDhbar (dn , dn) , expand , ' x ' ) ; 
grmap (DDhbar (dn , dn) , subs , EpsEqn , ' x ' ) ; 

grmap (DDhbar (dn, dn) , collect, 'yi' , ord, distributed, factor); 
grmap (DDhbar (dn,dn) , subs. Traces, 'xO; 

grmap (DDhbar (dn,dn) , collect, 'x^ , ord, distributed, factor); 
grdisplay (DDhbar (dn,dn) ) ; 

grdef ('Rhbar{a b} := 2*R{a "c b '•d}*hbar{c d}'); 
grcalc (Rhbar (dn , dn) ) ; 

grmap (Rhbar (dn,dn) , series, 'x' , e, 4); 
grmap (Rhbar (dn,dn) , convert, 'x', polynom); 
grmap (Rhbar (dn , dn) , expand , ' x ' ) ; 
grmap (Rhbar (dn, dn) , subs, EpsEqn, 'xO; 

grmap (Rhbar (dn, dn) , collect, 'x', ord, distributed, factor); 
grmap (Rhbar (dn, dn) , subs. Traces, 'x'); 

grmap (Rhbar (dn, dn) , collect, 'x', ord, distributed, factor); 
grdisplay (Rhbar (dn , dn) ) ; 

grdef ('LHS{(a b)} := DDhbar{a b} + Rhbar{a b}'); 
grcalc(LHS(dn,dn)) ; 

grmap (LHS (dn, dn) , series, 'x', e, 4); 
grmap (LHS (dn, dn) , convert, 'x' , polynom); 
grmap (LHS (dn , dn) , expand , ' x ' ) ; 
grmap ( LHS ( dn, dn ) , subs, EpsEqn, 'xO; 

grmap ( LHS ( dn, dn ) , collect, 'x' , ord, distributed, factor); 
grmap (LHS (dn , dn) , subs. Traces, 'x'); 

grmap (LHS (dn, dn) , collect, 'x' , ord, distributed, factor); 
grdisplay (LHS (dn , dn) ) ; 

grmap ( LHS ( dn, dn ) , subs, {e=l}, 'x'); 

grmap ( LHS ( dn, dn ) , collect, 'x\ ord, distributed, factor); 
grdisplay (LHS (dn , dn) ) ; 

grmap ( LHS ( dn, dn ) , collect, 'x\ [t, x, y, z, Exx, Eyy] , 
distributed, factor); 
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# Solution Of The Resulting Algebraic Equations Follows # 



#Final Verification Of The Solution # 

grdef('hl{(a b)} := kdelta{"$t a}*kdelta{'"$t b}*(-2)*mu*etadn{i 1}* 
AngM{''l}*BB{j k}*x{'-i}*x{^j}*x{'^k}/((x~2+y'-2+z'^2)''(3/2)) 
+ kdelta{^$t a}*kdelta{~$x b}*mu* 

(eps{$x i j}*etadn{l n}*AngM{"n}*etaup{"i "in}*EE{in k} 
*x{~j}*x{~l}*x{"k} 

- eps{$x i j}*AngM{''i}*EE{l k}*x{''j}*x{''l}*x{''k} 

+ 3*eps{i j l}*AngM{'^i}*EE{$x n}*etaup{''n ''j}*x{''l}*(x"2+y'^2+z'^2) ) 
/((x'-2+y'-2+z'-2)-(3/2)) 
+ kdelta-["$t a}*kdelta{"$y b}*mu* 

(eps{$y i j}*etadn{l n}*AngM{"n}*etaup{"i "m}*EE{m k} 
*x{'' j }*x{"l}*x{"k} 

- eps{$y i j}*AngM{''i}*EE{l k}*x{''j}*x{''l}*x{'"k} 

+ 3*eps{i j l}*AngM{~i}*EE{$y n}*etaup{''n '■j}*x{''l}*(x'^2+y"2+z"2) ) 
/((x-2+y"2+z-2)-(3/2)) 
+ kdelta{"$t a}*kdelta{~$z b}*mu* 

(eps{$z i j}*etadn{l n}*AngM{~n}*etaup{''i '■in}*EE{in k} 
*x{" j }*x{"l}*x{~k} 

- eps{$z i j}*AngM{"i}*EE{l k}*x-["j>*x-["l}*x{"k} 

+ 3*eps{i j l}*AngM{"i}*EE{$z n}*etaup{"n "j}*x{"l}*(x~2+y~2+z~2) ) 
/((x"2+y-2+z-2)-(3/2)) 
+ kdelta{''$x a}*kdelta{''$x b}*mu* 

(2*etadn{i k}*AngM{''k}*BB{$x $x}*x{''i}*(x'"2+y''2+z''2) 

+ (2/3)*AngM-["i}*BB{i j>*etadn-[$x k}*x{~k>*etadn{$x l}*x{"l}*x{'^j} 

+ 2*AngM-["i}*BB{$x i}*etadn-[$x k}*x{~k}* (x"2+y"2+z"2) 

- (2/3)*etadn{i k}*AngM{~k}*BB{$x j}*etadn{$x l}*x{"l}*x{~i}*x{~ j} 

- (10/3)*etadn{$x k}*AngM{~k}*BB{$x j}*x{''j}*(x''2+y'^2+z"2) 

- (2/3)*etadn{$x k}*AngM{''k}*BB{i j}*etadn{$x l}*x{''l}*x{"i}*x{'^j} 

- (2/3)*AngM-["i}*BB{i j}*x-["j}*(x"2+y"2+z"2) 

+ (2/3)*etadn{i l}*AngM{"l}*BB{j k}*x{~i}*x{" j}*x{"k}) 
((x-2+y~2+z~2)"(3/2)) 
+ kdelta-["$y a}*kdelta{~$y b}*mu* 

(2*etadn{i k}*AngM{''k}*BB{$y $y}*x{''i}*(x''2+y''2+z''2) 

+ (2/3)*AngM-["i}*BB{i j}*etadn-[$y k}*x{~k>*etadn{$y l}*x{"l}*x{'^j} 

+ 2*AngM-["i}*BB{$y i}*etadn-[$y k}*x{~k}* (x"2+y"2+z"2) 

- (2/3)*etadn{i k}*AngM{"k}*BB{$y j}*etadn{$y l}*x{~l}*x{~i}*x{~ j} 



98 



- (10/3)*etadn{$y k}*AngM{~k}*BB{$y j}*x{'' j}* (x"2+y~2+z"2) 

- (2/3)*etadn{$y k>*AngM{"k}*BB{i j}*etadn{$y l}*x{"l}*x{"i}*x{" j} 

- (2/3)*AngM-[~i}*BB{i j>*x{"j}*(x"2+y"2+z"2) 

+ (2/3)*etadn{i l}*AngM{"l}*BB{j k}*x{"i}*x{" j}*x{~k}) 
((x^2+y^2+z^2)~(3/2)) 
+ kdelta{''$z a}*kdelta{"$z b}*mu* 

(2*etadn{i k}*AngM{"k}*BB{$z $z}*x-["i}* (x"2+y"2+z"2) 

+ (2/3)*AngM-["i}*BB{i j}*etadn-[$z k}*x{"k}*etadn{$z l}*x{~l}*x{~ j} 

+ 2*AngM{"i}*BB{$z i}*etadn-[$z k}*x{~k}*(x"2+y"2+z"2) 

- (2/3)*etadn{i k}*AngM{''k}*BB{$z j}*etadn{$z l}*x{'"l}*x{'^i}*x{" j} 

- (10/3)*etadn{$z k}*AngM{''k}*BB{$z j}*x{''j}*(x'"2+y~2+z'^2) 

- (2/3)*etadn{$z k}*AngM{^k>*BB{i j}*etadn{$z l}*x{''l}*x{'^i}*x{"j} 

- (2/3)*AngM-["i}*BB{i j}*x{"j}*(x"2+y"2+z"2) 

+ (2/3)*etadn{i l}*AngM{"l}*BB{j k}*x{"i}*x{" j}*x{"k}) 
/((x^2+y^2+z~2)"(3/2)) 
+ kdelta{~$x a}*kdelta{'^$y b}*mu* 

(2*etadn{i k}*AngM{"k}*BB{$x $y}*x{"i}* (x"2+y"2+z"2) 

+ (2/3)*AngM-["i}*BB{i j}*etadn-[$x k}*x{"k}*etadn{$y l}*x{~l}*x{~ j} 

+ AngM-["i}*BB{i $x}*etadn{$y k}*x{~k}* (x~2+y~2+z~2) 

+ AngM{"i}*BB{i $y}*etadn{$x k}*x{'^k}*(x'^2+y"2+z"2) 

- (l/3)*etadn{i k}*AngM{'"k}*BB{j $x}*etadn{$y l}*x{'"l}*x{'^i}*x{"j} 

- (l/3)*etadn{i k}*AngM{"k}*BB-[j $y}*etadn{$x l}*x{"l}*x{~i}*x{"j} 

- (5/3)*AngM-[~k}*etadn{k $x}*BB-[$y j}*x{"j}*(x"2+y~2+z~2) 

- (5/3)*AngM{"k}*etadn{k $y}*BB{$x j}*x{" j}* (x"2+y"2+z~2) 

- (l/3)*AngM{'"k}*etadn{k $x}*etadn{$y l}*BB{i j}*x{''l}*x{"i}*x{'^j} 

- (l/3)*AngM{'"k}*etadn{k $y}*etadn{$x l}*BB{i j}*x{''l}*x{"i}*x{'^j}) 
/((x-2+y-2+z-2)-(3/2)) 

+ kdelta-["$x a}*kdelta{"$z b}*mu* 

(2*etadn{i k}*AngM{"k}*BB{$x $z}*x{"i}* (x"2+y"2+z"2) 

+ (2/3)*AngM{'"i}*BB{i j}*etadii{$x k}*x{'^k}*etadn{$z l}*x{'^l}*x{'^j} 

+ AngM{''i}*BB{i $x}*etadn{$z k}*x{'^k}*(x"2+y~2+z~2) 

+ AngM-["i}*BB{i $z}*etadn{$x k}*x{~k}* (x~2+y~2+z~2) 

- (l/3)*etadn{i k}*AngM{~k}*BB-[j $x}*etadn{$z l}*x{"l}*x{~i}*x{" j} 

- (l/3)*etadn{i k}*AngM{"k}*BB{j $z}*etadn{$x l}*x{"l}*x{~i}*x{~ j} 

- (5/3)*AngM{''k}*etadn{k $x}*BB{$z j}*x{''j}*(x''2+y''2+z"2) 

- (5/3)*AngM{''k}*etadn{k $z}*BB{$x j}*x{''j}*(x''2+y''2+z"2) 

- (l/3)*AngM-["k}*etadn{k $x}*etadn{$z l}*BB{i j}*x{"l}*x{~i}*x{~ j} 

- (l/3)*AngM{"k}*etadn{k $z}*etadn{$x l}*BB{i j}*x{"l}*x{~i}*x{~ j}) 
/((x'-2+y-2+z-2)-(3/2)) 

+ kdelta-["$y a}*kdelta{~$z b}*mu* 

(2*etadn{i k}*AngM{''k}*BB{$y $z}*x{''i}*(x''2+y''2+z''2) 

+ (2/3)*AngM-["i}*BB{i j}*etadn-[$y k}*x{~k>*etadn{$z l}*x{"l}*x{'^j} 

+ AngM-["i}*BB{i $y}*etadn{$z k}*x{~k}* (x"2+y"2+z"2) 

+ AngM{"i}*BB{i $z}*etadn{$y k}*x{~k}* (x~2+y~2+z~2) 
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- (l/3)*etadn{i k}*AngM{^k}*BB{j $y>*etadn{$z l}*x{^l}*x{~i}*x{" j} 

- (l/3)*etadn{i k}*AngM{"k}*BB{j $z}*etadn{$y l}*x{"l}*x{~i}*x{" j} 

- (5/3)*AngM-["k}*etadn{k $y}*BB-[$z j}*x{"j}*(x~2+y~2+z~2) 

- (5/3)*AngM{"k}*etadn{k $z}*BB-[$y j}*x{" j}* (x~2+y~2+z~2) 

- (l/3)*AngM{'"k}*etadn{k $y}*etadn{$z l}*BB{i j}*x{''l}*x{"i}*x{'^j} 

- (l/3)*AngM{'"k}*etadn{k $z}*etadn{$y l}*BB{i j}*x{''l}*x{"i}*x{'^j}) 
/((x'-2+y'-2+z-2)-(3/2)) ' ): 

grcalc (hi (dn , dn) ) : 

grdef ('hh{(a b)} := hO{a b} + hl{a b}'): 
grcalc (hh (dn , dn) ) : 

grdef ( 'htrace := hh{~c c}'): 

grcalc (htrace) ; grmap (htrace, series, 'x' , e, 4): 
grmap (htrace , convert, 'x' , polynom) : 
grmap (htrace , expand, 'x'): 
grmap (htrace , subs, EpsEqn, 'x'): 

grmap(htrace , collect, 'x', ord, distributed, factor): 
grmap (htrace , subs. Traces, 'x'): 

grmap (htrace , collect, 'x' , ord, distributed, factor): 
grdisplay (htrace) : 

grdef ('hbar-[(a b)} := hh{a b} - (l/2)*g{a b}*htrace'): 
grcalc (hbar(dn,dn)) : 

grmap (hbar (dn, dn) , series, 'x' , e, 4): 
grmap (hbar (dn, dn) , convert, 'x' , polynom): 
grmap (hbar (dn , dn) , expand , ' x ' ) : 
grmap (hbar (dn, dn) , subs, EpsEqn, 'x'): 

grmap(hbar (dn,dn) , collect, 'x' , ord, distributed, factor): 
grmap (hbar (dn, dn) , subs. Traces, 'x'): 

grmap (hbar (dn, dn) , collect, 'x' , ord, distributed, factor): 
grdisplay (hbar (dn, dn) ) : 

grcalc (hbar (dn , dn , cdn) ) : 

grmap (hbar (dn, dn, cdn) , series, 'x' , e, 4): 

grmap (hbar (dn,dn, cdn) , convert, 'x', polynom): 
grmap (hbar (dn , dn , cdn) , expand , ' x ' ) : 

grmap(hbar (dn,dn, cdn) , subs, EpsEqn, 'x'): 

grmap (hbar (dn,dn, cdn) , collect, 'x', ord, distributed, factor): 

grmap(hbar (dn,dn, cdn) , subs. Traces, 'x'): 

grmap (hbar (dn,dn, cdn) , collect, 'x', ord, distributed, factor): 
grdisplay (hbar (dn , dn , cdn) ) : 

grdef ('Dhbar{(a b) "c} := hbar{a b Tc}'): 
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grcalc (Dhbar (dn , dn , up) ) : 

grmap(Dhbar(dn,dn,up) , series, 'x' , e, 4): 
grmap (Dhbar (dn,dn, up) , convert, 'x', polynom) : 
grmap (Dhbar (dn , dn , up) , expand , ' x ' ) : 
grmap (Dhbar (dn,dn, up) , subs, EpsEqn, 'x'): 

grmap (Dhbar (dn,dn, up) , collect, 'x', ord, distributed, factor): 

grmap (Dhbar (dn,dn, up) , subs. Traces, 'x'): 

grmap (Dhbar (dn,dn, up) , collect, 'x', ord, distributed, factor): 
grdisplay (Dhbar (dn , dn , up) ) : 

grdef ('DDhbar{(a b)} := Dhbar{a b ~c ;c}'): 
grcalc (DDhbar(dn,dn)) : 

grmap (DDhbar (dn,dn) , series, 'x', e, 4): 
grmap (DDhbar (dn,dn) , convert, 'x' , polynom): 
grmap (DDhbar (dn , dn) , expand , ' x ' ) : 
grmap (DDhbar (dn , dn) , subs , EpsEqn , ' x ' ) : 

grmap ( DDhb ar (dn,dn) , collect, 'x' , ord, distributed, factor): 
grmap ( DDhb ar (dn,dn) , subs. Traces, 'x'): 

grmap (DDhbar (dn,dn) , collect, 'x' , ord, distributed, factor): 
grdisplay (DDhbar (dn, dn) ) : 

grdef ('Rhbar{a b} := 2*R{a '"c b '•d}*hbar{c d}'): 
grcalc (Rhbar(dn,dn)) : 

grmap (Rhbar (dn, dn) , series, 'x', e, 4): 
grmap (Rhbar (dn, dn) , convert, 'x', polynom): 
grmap (Rhbar (dn , dn) , expand , ' x ' ) : 
grmap (Rhbar (dn,dn) , subs, EpsEqn, 'x'): 

grmap (Rhbar (dn, dn) , collect, 'x', ord, distributed, factor): 
grmap (Rhbar (dn, dn) , subs. Traces, 'x'): 

grmap (Rhbar (dn,dn) , collect, 'x', ord, distributed, factor): 
grdisplay (Rhbar (dn , dn) ) ; 

grdef ('LHS{(a b)} := DDhbar{a b} + Rhbar{a b}'): 
grcalc (LHS(dn,dn)) : 

grmap (LHS (dn , dn) , series, 'x', e, 4): 
grmap (LHS (dn, dn) , convert, 'x' , polynom): 
grmap (LHS (dn , dn) , expand , ' x ' ) : 
grmap ( LHS ( dn, dn ) , subs, EpsEqn, 'x'): 

grmap ( LHS ( dn, dn ) , collect, 'x' , ord, distributed, factor): 
grmap ( LHS ( dn, dn ) , subs. Traces, 'x'): 

grmap (LHS (dn, dn) , collect, 'x' , ord, distributed, factor): 

grdisplay (LHS (dn,dn)) ; 

grmap ( LHS ( dn, dn ) , subs, {e=l}, 'x'): 

grmap (LHS (dn, dn) , collect, 'x' , ordtxyz, distributed, factor): 



grdisplay (LHS (dn , dn) ) 



APPENDIX B 
DECOMPOSITION OF p^x'^ 



A detailed description of tiie spiierical iiarmonic decomposition of the product 
^X~'^ is presented in this appendix, where p is any integer and q can be any real 
number, even though in every case in which this decomposition will be used q will be 
either an integer or a half-integer. The notation used is 

p2 ^ ''o^' +2r' /''~l^ v{l-cosQ) (B.l) 
^ To - 2M ° To - 3M ^ ^ ^ ^ ^ 

X = 1 ^^^TTsin^^ (B.2) 

^ ro-2M ^ ^ 

A = r-To. (B.3) 
In order to simplify the equations, the following symbols are used: 

^ = ^ (B-4) 

ro-2M ^ ^ 

° r„ - 3M ^ ^ 

^ A-(r.-3M) 

2r„(r„-2M)2 ^ ^ 

^2 ^ A2(r,-3M) 7^ 
2ro(r„ - 2M)2x x ' 

Using Equations (B.5)-(B.8) one gets 

2 

p = B^x^ (— + l-cose)^ =5^X^ (5^ + 1-cose)^ (B.9) 
X 

pf;^^5 ^ x^^"^ {— + 'i--cose)^ ^ B2 x^"^ {5'^ + l-cose)^. (B.IO) 

X 



(B.8) 
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For the decomposition then it is assumed that 

= x^"^ (— + 1 - cos e) ^ 

oo / (B.ll) 

1=0 m=-l 

and the coefficients Ef'^ are calculated by multiplying both sides of the above equation 
with the complex-conjugate spherical harmonics y5^(0, $) and integrating over = 
[0, tt] and $ = [0, 2n]. Using the results given in the discussion of spherical harmonics 
in [38], the following expression for the coefficients is obtained 



(B.12) 

where u = cos O and P"\u) are the Associated Legendrc polynomials. In the following 
it is also assumed that Pl^^^{u) = Pi{u), which are the Legendre polynomials. It 
should be noted that, because only even powers of sin $ show up in the expression for 
the coefficients Ef'^, only the even values of m contribute to the sum (B.ll) and for 
all odd m the coefficients Ef'^ are equal to 0. In all the equations that follow, m is 
assumed to be even. 

First, the integral 

IL = £ du{S' + 1 - uf^Pr\u) (B.13) 

is calculated. Assuming that m is positive. Equation 8.911 of [42] gives the m- 
derivative of the Legendre polynomials 

d^'Piju) ^ 1 ^ (-l)^r(2/-2A; + l) 

where r(a;) is the usual F-function and [|] stands for the largest integer smaller or 
equal to |. Then, Equation 8.810 of [42] gives for the associated Legendre polynomials, 
for positive m 



p^r,,! - n _,,2^f 1 V (-1) F(2/-2/c + l) l-2k-m /-oip-N 
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where the right-hand side is a polynomial, since m is even. The case of negative m 
can be dealt with if one takes into account the fact that for any m: 

^r"W = ^|^/'r'M (B.ia) 

which can be easily derived from Equation 8.752 of [42]. Finally, the expression of the 
associated Legendre polynomials that is used in the calculation of the integral If^ is 
found to be 

p-m. x _ r(/-|m| + l) 1 ^ 

r(Z + m + l) 2^'' 

(-l)^r(2/-2/c + l) ^,_,,_|^| 



^ r(A; + l)r(/ -k + l)r(/ - 2A; - |m| + 1) 



(B.17) 



which is valid for all even m, positive or negative. Substituting Equation (B.17) into 
Equation (B.13) gives 



_ r(/-H + i) i^ (-im2/-2fc + i) 

T{l + m + 1) 2' f^^T{k + l)r(/ - A; + l)r(/ -2k- \m\ + 1) 



The integral of Equation (B.18) can be calculated if the term (1 — u'^) 2 is expanded 
in powers of u, specifically 



M , , (B.19) 

t^r(^ + i)r(^-^ + 1)7-1 

The task, then, shifts to the calculation of the integral 

I{a,b)^ J duu%5^ + l-u)'^ (B.20) 
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where a and b are integers, with a > 0. Performing integration by parts once results 
in 



'(a,6) 



-1 
1 



^ + 1 

2 



is 



+ 1 



2\^ + l 



;-ir^(2 + 5^)^^+ 2« 



- + 1 

2 ~ 



b + 2 



- (52)1+1 + (_l)a(2 + 5^)1+1 



b + 2 

2a ^ 



b + 2 



'(a-1,6+2) 



which is a recursion relation for /(„ 5) . Using the fact that 

2 



'(0,6) 



b + 2 



{5')-^+^ + {2 + 6')-^- 



(B.21) 



(B.22) 



it can be proven by induction that the general integral I(^a,b) is equal to 

2"+i r(a + 1) 



^^"'"^ S + + 4) . . . [6 + 2(n + 1)] r{a-n + 1) 



-l)»-"(2 + 52)§+'*+i 



;B.23) 



Finally, using Equations (B.23) and (B.19) into Equation (B.18), the integral If!^ ends 
up being 



— 



r{i- \m\ + l) 1 



E 



;-i)'^r(2/-2fc + i) 



r(l + m + 1) 2^ ^ T{k + 1)T{1 -k + 1)T{1 -2k- \m\ + 1) 



X 



Jm| 
2 

E 



r(M + i) 



Z-|m|-2fc+2i/ 



E 



)n+l 



^r(i/ + l)r(^-i/ + l) ^ (p + 2)(p + 4)...[p + 2(n + l)] 

[(-l)'-"(2 + (52)f+«+i _ (<52)f+«+i 



X 



(B.24) 



r(/- 




m 


-2k + 2u + l) 


r(/- 


m 




-2k + 2u-n + l) 



A short comment about the term (52)f+"+^ should be made at this point. If p is 
an odd integer, that term implies that the square root of S"^, or of 72, must be taken. 
One thing that may be ambiguous is the sign of the square root, + or — , that must 
be considered. That point is discussed in Chapter 6 of this dissertation, where the 
physical interpretation of the quantities becomes clear and the correct sign to use can 
be easily recognized. For the calculations in this appendix it is enough to interpret 
the term (5^)2 as the appropriate square root of 5'^. 
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Continuing the calculation of the coefficients Ef'^, substituting Equation (B.24) 
into Equation (B.12) shows that the next two integrals to be done are 



I, = (-1)^-" / d^e-'^^^x'^'^C^ + 
Jo 

= (-1)'-" [ 
Jo 

Jo 



X 



:B.25) 



p2n 



Jo 



(B.26) 



X 

and their difference must be multiplied by the appropriate factor given by Equation 
(B.24). For both these integrals, the exponential e~*"** must be written in terms of 
the variable x- general 

g-im* ^ cos(m$) - i sin(m<l') (B.27) 

but, since only even m's are considered and the remaining of the integrand is a function 
of sin^ $, the sin(m$) part of the exponential gives once the integration is performed. 
So, as far as the integrals are concerned, the exponential can be replaced by the 
cos(m$). Then Equation 1.331 of [42] gives that for any even m, positive or negative, 

cos(m$) = y ^„ r cosl"^l-2^ $ sin^^ $. (B.28) 

^ ' 4^r(|m| -2A + l)r(2A + l) ^ ^ 

A — 

The differential (i$ is written in terms of the differential dx 

= ^ . (B.29) 

2C sin $ cos $ 



Also 



sin2$ = 1—^ (B.30) 



cos^$ = l-^—^ (B.31) 

and caution should be used when the square root of these expressions is taken. Specif- 
ically, for $ in [0, |) or in [tt, ^), sin $ cos $ > 0, so the equation 

sin $ cos $ = ^(^)(l-^^^ (B.32) 
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should be used. On the other hand, for $ in [f,7r) or in [^,27r), sin $ cos $ < 0, so 
the equation 

sin $ cos $ = - a/(^)(i-^^ (B.33) 



should be used. 

In 7i, the factor (2 + 2_ jg first written as an infinite sum 



I ^%g+n+l _ r(^+n+l) 2+„_^+i7^^ f-noA] 



The integral becomes 



I m| 



T -( r(H + i)(-i)^ Mf. r(f + n + i) 

' ^ ^ ^r(|m|-2A + l)r(2A + l) ^T{s + l)T{^ + n-s + l) 



A=0 



Jl-C 

(B.35) 

Changing the variable of integration from % to y = and using the definition of 
the hypergeometric functions 

where: 3^(7) > 3?(/9) > 0, result in the final expression for the integral 



I m| 



X 



= 0F(-i)'-" V(-i)^ y 7'^2f - 

r(A + l)r(M_A + l)^^ (B_37) 



r(| + n + l) ^ n ^1 |m| ^ ^, 

r^(-^ + ? + s>'^+ n; V + l;C')■ 



^(s+l)^(|+r^-s+l) ' 2 ^ ' 2' 2 

To calculate the second integral I2 one follows the same procedure that was used 
for Ii, the only difference being that the step of expanding the factor (2 + ^^^f+i^+i 
is now not necessary. The result is 

|m| 

2 p/ \m\ + l \ ^ I I 

^ r(A + l)r(H_A + l)^' ^ ' 2' 2 

(B.38) 
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The final expression for the coefficients Ef'^{'y^, C) can be obtained by combining 
Equations (B.24), (B.37), (B.38) and (B.12). It is 



E»^,^c) = JHi±i ni-\m\ + i) — 

• V IT [r(/-m+l)r(( + m + l)]i 

|m| 

^r(A + i)r(H_A + i)'' 

E/_i \k2i-2k r{i — k + ^) 
fc=o r(A; + l)r(Z-2A;- |m| + 1) 

|m| 

r(H + i) 

^r(^ + i)r(M_^ + i) 

/-|m|-2A;+2!/ 

E 

n=0 

_2(f)f+«+iF(g + n + l,A + ^;^ + l;C) + 

c_i y+" ^,,2s r(| + n + i) E+^_^+2 
^ ^ r(s + i)r(| + n-s + i) 



2" r(/ - 




m 


- 2A; + 2i/ + 1) 


(p + 2)(p + 4)...[p + 2(n + l)]r(Z-l 


m 




-2k + 2u-n + l) 



X 



F(-| + g + .,A+i;^ + l;C) 



;B.39) 



It is reminded that this expression is vahd for all even m's, and that for odd m's these 
coefficients are equal to zero. 

If this expression for the coeffients is used, the equations become very long 
and their dependence on 7^, which is of particular interest for the calculation of the 
regularization parameters, becomes unclear. In order for the expressions that appear 
in this dissertation to look simpler, the following notation is used for the three finite 
sums over k, u and A which appear in the previous equation and which have no 
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7^-dependence 



2/ + 1 r(/-|m| + l) 

TT [r(Z-m+l)r(/ + m + l)]5 



^ V{k + l)r(Z - 2A; - |m| + 1) ^ r(i/ + i)r(^ - 1/ + 1) 



fe=0 

|m| 

V(-i)^ LLjzi 

r(A + i)r(M_A + i)- 

(B.40) 



It is implied that {Ylm v\) ^ multiplicative factor and any quantities appearing to 

(a) 

2' 2 



the right of it can depend on /c, v or A. Also, the notation F^^^ is used for the 

hypergeometric function F(a,A+|;^ + l;C). 
Finally 



OT^c)=(E) E 



Z-|m|-2fe+2i/ 



,,,,, n=o (p + 2)(p + 4)...[p + 2(n+l)] 



X 



r(/- 




m 


- 2A; + 2i/ + 1) 


ni-\ 


\m 







E 2g ( + n + 1) g^^_^^2 c^A.m 

^ r(s + i)r(f + n-S + l) ^(-1+9+-) f 



s=0 



APPENDIX C 
NUMERICAL CODE FOR THE RETARDED FIELD 



#include <stdlib.h> 
#include <iostreani.h> 
#include <iomanip.h> 
#include <f stream. h> 
#include <assert.h> 
#include <string.h> 
#include <stdio.h> 
#include <complex.h> 
#include <matli.h> 



// Convenient macros for debugging 

#define SHOW(a) " "«#a«" = "«a«" " 

#define PRT(a) cerr « " "«#a«" = " « a « endl; 

#define PRINT(a) cerr «__LINE__« " 1" «__FILE__«": "« 

#a «" = "«a« endl; 
#def ine WHEREAMI __LINE__ « " I " « __FILE__ « " : " 
#define TRACE cerr « __LINE__ « " 1" « __FILE__ « ": " 
#define TRACK cerr « __LINE__ « " 1" « __FILE__ « endl; 



typedef complex<double> cmplx ; 



int kmaXjkount; 
double *xp,dxsav; 
cmplx **yp; 



cmplx muCint 1, int m) ; 
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Ill 



void f ksehCcmplx *fl, cmplx *f2, cmplx *f3, int k) ; 
void fksinf (cmplx *fl, cmplx *f2, cmplx *f3, int k) ; 
cmplx coefkeh(int k) ; 
cmplx coefkinf(int k) ; 

void rkqsCcmplx y[], cmplx dydr[], int n, double *r, double htry, 
double eps, double yscal[], double *hdid, double *hnext, 
void (*derivs) (double, cmplx [], cmplx [])) ; 

void rkckCcmplx y[], cmplx dydr[], int n, double r, double h, 
cmplx yout [] , cmplx yerr[], void (*derivs) (double, 
cmplx [] , cmplx [] ) ) ; 

void derivs (double r, cmplx y[], cmplx dydr[]) ; 

void factor (cmplxfe coefr, cmplx& coefq, double r) ; 

void strline (double xl, double x2, cmplx yl, cmplx y2, 
double R,cmplx& yR) ; 

void odeint(cmplx ystart[],int nvar, double xl, double x2, double eps, 
double hi, double hmin, int *nok, int *nbad, 
cmplx& ylastl, cmplxfe ylast2, 
void (*derivs) (double, cmplx [], cmplx [] ) , 
void (*rkqs) (cmplx [], cmplx [] , int, double *, 
double, double, double [] , double *, double *, 
void (*) (double, cmplx [], cmplx []))) ; 



const cmplx 1(0,1); 

namespace ScalarWave{ 
double M = 1 ; 
int 1; 
int m; 

double omega; 
double w; 

double R = 20*M; 

double dtdtau = 1/sqrt (1-(3*M/R) ) ; 

} 



// Mass of the Black Hole 

// Spherical Harmonic Index 

// Spherical Harmonic Index 

// Frequency of Radiation 

// Angular Frequency of Orbit 

// Radius of the Orbit 

// t-Component of 4-velocity 
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#include "rkqs.c" 
#include "rkck.c" 
#include "nrutil.c" 
#include "odeintl.c" 
#include "sum.c" 
#include "nrutil.h" 
#include " sphharm . c " 
#include "zero.cpp" 



int inain(){ 

using namespace ScalarWave; 
M = 1.0 ; 

w = sqrt(M/(R*R*R)) ; 
double pi = 4.0*atan(1.0) ; 

char* buff; 

buff = new char [513]; 

char bh; 

int nvar =2; // Number of First Order Differential Equations, 
int nok = 0; 
int nbad = 0; 
int Imax = 31 ; 

cmplx ystart [3] ; 
double mulm; 

cmplx yLRl, yLR2, yRRl, yRR2; 
cmplx Wl, W2, W3, W4 ; 

cout « endl; 

cout « "Radius of the orbit:" « SHOW(R) « endl; 
of stream out ( "DPsi20 . dat " ) ; 
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out « "1 m mulm yLN2 yRN2" « endl; 



// Initialization of the Self -Force Contributions. 

double Flmt =0; 
double Fit =0; 
double Ft =0; 
double FlmrL =0; 
double FlrL =0; 
double FrL =0; 
double FlmrR =0; 
double FlrR =0; 
double FrR =0; 



for (1=0; Klmax; 1++) { 



for(m=l ; in>-l ; in=in-2) { 



if (m!=0) { 
omega = in*w ; 



cout << endl ; 

cout « SHOW(l) « SHOW(m) « endl; 

cmplx bl = (1*1+1+1) / (2*M+8*I*omega*M*M) ; 
cmplx al = -I*l*(l+l)/(2*omega) ; 

double Reh = l/(100*abs(bl)) + 2*M ; 

// Starting Radius Close to Event Horizon 
double Rinf = 100*abs(al) ; 

// Starting Radius Close to Infinity 

// Integrate in the Region (Rs, R] to Get yL(r) . 



double h = 1 . Oe-5 ; 
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double hmin = l.Oe-10 ; 
double eps = l.Oe-12 ; 

kmax = 0; 
dxsav = 0.0; 

// Boundary Conditions at the Event Horizon : 

cmplx expeh = 1 ; 
cmplx Dexpeh = 0; 
coefkeh(O) ; 
double Xn = 1; 
double last = l.e+33; 



forCint k=l; k<100; k++) { 

Xn *= (Reh-2*M) ; 
cmplx bk = coefkeh(k); 
cmplx next = bk*Xn; 

if (expeh + next == expeh) { 

cout « "Machine precision in event -horizon loop" 

« endl; 

cout « setprecision(13) « SHOW(k) « SHOW(expeh) 

« SHOW (next) « endl; 
breai ; 

} 



if (abs(next) > last) { 

cerr « "No convergence in event -horizon loop" 
« endl; 

cerr « setprecision(13) « SHOW(k) « SHOW(expeh) 

« SHOW (next) « endl; 
break; 



} 
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last = abs(next); 

expeh += next ; 

Dexpeh += k*next/ (Reh-2*M) ; 

if (k < 5) { 

cout « setprecision(13) « SHOW(k) « SHOW(bk) 
« SHOW(expeh) « endl; 

} 

} 



expeh =expeh* 

expCcmplxCO.O, omega*Reh + omega*2.0*M*log(Reh/(2*M)-l) )); 
Dexpeh = exp(cmplx(0.0, omega*Reh + 

omega*2 . 0*M*log (Reh/ (2*M) -1) ) ) *Dexpeh 
+ expeh*cmplx(0 , omega) *Reh/ (Reh-2*M) ; 

ystart [1] =expeh/Reh; 

ystart [2] =Dexpeh/Reh - expeh/ (Reh*Reh) ; 

cout « setprecision(13) « SHOW(ystart [1] ) 
« SHOW (ystart [2]) « endl; 

// Wronskian at the event horizon. 

Wl = Reh* (Reh-2*M) * ( ystart [1] *conj (ystart [2] ) 

- ystart [2] *conj (ystart [1] ) ) ; 

cout « SHOW(Wl) « SH0W(-2*I*oinega) « endl ; 

// Integration of the differential equation for Psi. 

odeint (ystart , nvar , Reh , R , eps , h , hmin , &nok , fenbad , yLRl , yLR2 , 
(*derivs) , (*rkqs)) ; 



116 



// Wronskian at the orbit . 

W2 = R*(R-2*M)* (yLRl*conj (yLR2) - yLR2*conj (yLRl) ) ; 
//Solution before the matching. 

cout « setprecision(13) « SHOW(Reh) « SHOW(yLRl) 

« SH0W(yLR2) « endl; 
cout « setprecision(13) « SHOW(Wl) « SH0W(W2) 

« SH0W(-2*I*omega) « endl; 

//Integrate in the region [R.Rinf ) to Get yR(r) . 

h = -l.Oe-5; 
hmin = l.e-10; 
eps = l.e-12; 
kmax = 0; 
dxsav = 0.0; 



//Boundary Conditions at Infinity: 

cmplx expinf = 1; 
cmplx Dexpinf = 0; 
coef kinf (0) ; 
Xn = 1; 

last = l.e+33 ; 



forCint k=l; k<100; k++) { 
Xn *= Rinf; 

cmplx ak = coef kinf (k) ; 
cmplx next = ak/Xn; 

if (expinf + next ==expinf) { 

cout « "Machine precision in infinity loop" « endl; 
cout « setprecision(13) « SHOW(k) « SHOW(expinf) 
« SHOW (next) « endl; 
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break; 

} 

if (abs(next) > last) { 

cout « "No convergence in infinity loop" << endl; 
cout « setprecision(13) « SHOW(k) « SHOW(expinf) 

« SHOW (next) « endl; 
break; 

} 



last = abs(next); 

expinf += next; 

Dexpinf += -k*next/Rinf ; 

if (k<5) { 

cout « setprecision(13) « SHOW(k) « SHOW(ak) 
« SHOW (expinf) « endl; 

} 



} 

double rstar = Rinf + 2*M*log(Rinf /(2*M)-1) ; 

expinf = exp(l*omega*rstar) *expinf ; 
Dexpinf = I*oinega*Rinf *expinf /(Rinf-2*M) 
+ exp(I*omega*rstar)*Dexpinf ; 

ystart[l] = expinf /Rinf; 

ystart [2] = Dexpinf /Rinf - expinf /(Rinf *Rinf) ; 

cout « setprecision(13) « SHOW (ystart [1] ) 
« SHOW (ystart [2]) « endl; 
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//Wronskian at infinity 

W3 = Rinf*(Rinf-2*M)*(ystart[l]*conj (ystart[2]) 

-ystart [2] *conj (ystart [1] ) ) ; 



//Integration of the differential equation for Psi. 

odeint (ystcirt , nvar , Rinf , R , eps , h , hmin , &nok , fenbad , yRRl , yRR2 , 
(*derivs) , (*rkqs) ) ; 



/ /Wronskian at the orbit . 

W4 = R* (R-2*M) * (yRRl*conj (yRR2) -yRR2*conj (yRRl) ) ; 



//Solution before the matching. 

cout « setprecision(13) « SHOW(Rinf) « SHOW(yRRl) 

« SH0W(yRR2) « endl; 
cout « setprecision(13) « SH0W(W3) « SH0W(W4) 

« SH0W(-2*I*omega)« endl ; 



//Matching of the solutions at the orbit. 

cmplx A = -mu(l,m) * yRRl /(R-2*M)/ (yLRl*yRR2-yRRl*yLR2) ; 
cmplx B = -inu(l,m) * yLRl /(R-2*M)/ (yLRl*yRR2-yRRl*yLR2) ; 



//Normalized solution. 

cmplx yLNl = yLRl* A ; 

cmplx yLN2 = yLR2*A ; 

cmplx yRNl = yRRl*B ; 

cmplx yRN2 = yRR2*B ; 



cout « SHOW(yLRl) « SH0W(yLR2) « endl; 
cout « SHOW (yRRl) « SH0W(yRR2) « endl; 
cout « SHOW (A) « SHOW(B) « endl; 
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cout « SHOW(yLNl) « SH0W(yLN2) « endl; 

cout « SHOW(yRNl) « SH0W(yRN2) « endl « endl; 



// Contribution to the self -force 

Flmt = 2*real(I*omega*yLNl*sharm(l,m,pi/2.0,0.0)) ; 
FlmrL = 2*real(yLN2*sharm(l,m,pi/2.0,0.0)) ; 
FlmrR = 2*real(yRN2*sharm(l,m,pi/2.0,0.0)) ; 



out « 1 «" "« m «" "« mu(l,in) «" "« yLN2 
«" "« yRN2 « endl; 



} 



if (m==0) { 

cmplx PsiH, dPsiH, Psil, dPsil; 



cout « endl; 

cout « SHOW(l) « SHOW(m) « endl; 



getZeroPsiH(&PsiH,&dPsiH) ; 
getZeroPsiI(&PsiI,&dPsiI) ; 



cmplx wronsk = R*(R-2.0*M)*(PsiH*dPsiI-PsiI*dPsiH) ; 



cmplx PsiHN = -mu(l,m) * PsiH*R*PsiI/wronsk; 
cmplx PsilN = -mu(l,m) * PsiI*R*PsiH/wronsk; 
cmplx DPsiH = -mu(l ,m) *dPsiH*R*PsiI/wronsk; 
cmplx DPsil = -mu(l,m)*dPsiI*R*PsiH/wronsk; 



Flmt = real (I*omega*PsiHN*sharm(l,m, pi/2. 0, 0. 0) ) ; 
FlmrL = real (DPsiH*sharm(l ,m, pi/2 . , . 0) ) ; 
FlmrR = real (DPsiI*sharm(l,m, pi/2. 0, 0. 0) ) ; 
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cout « SHOW(DPsiH) « SHOW(DPsil) « endl « endl; 
out "«m«" "«mu(l,m)«" "«DPsiH«" " 

«DPsiI« endl; 

} 



Fit += Flmt; 
FlrL += FlmrL; 
FlrR += FlmrR; 

} 



Ft += Fit; 
FrL += FlrL; 
FrR += FlrR; 



} 



out . close () ; 

cout « "File DPsi20.dat closed" « endl; 

} 



// Factor mu_(l,m). 

cmplx mu(int 1, int m){ 

using namespace ScalarWave; 

double mmu = 1.0; 

double pi = 4.0*atan(1.0) ; 



cout «SHOW (R) «SHOW (M) «SHOW (1 ) «SHOW (m) 
«SHOW ( sharm (l,m,pi/2,0))« endl ; 
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cmplx mulm = 4.0*pi*mmu * conj (sharm(l,m, pi/2. 0,0.0)) 
*sqrt(l-(3*M/R))/R; 

return mulm ; 



// Derivatives, For Use in rkqs.c . 

void derivs (double r, cmplx y[], cmplx dydr[]) { 

cmplx coefr,coefq ; 

f actorCcoef r , coefq, r) ; 

dydr[l] = y[2] ; // dy/dr 

dydr[2] = coefr*y[l] - coefq * y[2] ; // d"2y/dr'^2 

return ; 

} 



// coefr(r) and coefq(r) 

void factor (cmplxfe coefr, cmplxfe coefq, double r) { 

using namespace ScalarWave; 

int 12 = 1*(1+1); 
double r2 = r*r; 
double rm = r-2*M; 

coefr = -( omega*omega*r2/rm/rm - 12/r/(r-2.0*M)) ; 
coefq = 2.0*(r-M)/r/rm ; 



} 



// Spherical Harmonics. 



cmplx sharm(int 1, int m, double theta, double phi){ 

cmplx I = cmplx (0, 1 .0) ; 
double c = cos (theta) ; 

cmplx sh = exp(I*abs(m)*phi) * plgndr(l,abs(m) ,c) ; 
float fact = 

sqrt ( (2*1+1) *factorial (1-abs (m) ) 

/ (4 . 0* (4 . 0*atan ( 1.0))) /factorial (1+abs (m) ) ) ; 

sh = sh*fact; 

if(m<0) sh = pow(-l,abs(m))*conj (sh) ; 
return sh; 

} 



// Terms in the expansion of the intial guess for the solution 
// at the event horizon and at infinity 

void fksinf (cmplx *fl, cmplx *f2, cmplx *f3, int k){ 
using namespace ScalarWave; 

*fl = cmplx(-k*(k-l)+l*(l+l) , 4*omega*M*(k-l) ) 

/ cmplx (0.0, 2*omega*k) ; 
*f2 = cmplx( 0.0, -M*( (k-l)*(2*k-3) - 1*(1+1) ) ) /omega/k 
*f3 = cmplx ( 0.0, 2*M*M*(k-2)*(k-2) /omega/k ); 



void fkseh(cmplx *fl, cmplx *f2, cmplx *f3, int k){ 
using namespace ScalarWave; 
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*fl = cmplxC - (2*k-3)*(k-l)+ (1*1+1+1) , -12*omega*M*(k-l) ) 

/2.0/M/ cmplx(k*k, 4*omega*M*k) ; 
*f2 = cmplxC -(k-2)*(k-3)+l*(l+l) , -12*omega*M* (k-2) ) 

/4.0/M/M/ cmplx(k*k, 4*omega*M*k) ; 
*f3 = cmplxC 0.0, -omega* (k-3) ) 

/ 2.0/M/M/ cmplx(k*k, 4*omega*M*k) ; 

} 



cmplx coefkehCint KK) { 

using namespace ScalarWave; 

static cmplx bk_l = 1 ; 
static cmplx bk_2 = 0; 
static cmplx bk_3 = 0; 
static int k = 0; 

if (KK == 0) { 
bk_l = 1; 
bk_2 = 0; 
bk_3 = 0; 
k = 0; 

} 

if (KK != k) { 

cerr « "sum: : coefkeh called out of order, k = " 

« KK « " expected " « k « endl; 
assert (KK == k) ; 

} 

cmplx fl = 0; 
cmplx f2 = 0; 
cmplx f3 = 0; 
cmplx ret; 

if (k == 0) { 
k++; 

return 1; 
} else {. 

fkseh(&fl, &f2, &f3, k) ; 



ret = fl*bk_l + f2*bk_2 + f3*bk_3; 

} 



bk_3 = bk_2; 
bk_2 = bk_l; 
bk_l = ret; 

k++; 

return ret; 

} 



cmplx coefkinf(int KK) { 
using namespace ScalarWave; 



static 


cmplx 


ak_l 


= 1 


static 


cmplx 


ak_2 


= 


static 


cmplx 


ak_3 


= 


static 


int k 


= 0; 




if (KK 


== 0) 


{ 




ak_l 


= 1; 






ak_2 


= 0; 






ak_3 


= 0; 






k = 0; 






} 








if (KK 


!= k) 


i 





cerr « "sum: : coef kinf called out of order, k 

« KK « " expected " « k « endl; 
assert (KK == k) ; 

} 

cmplx fl, f2, f3, ret; 

if (k == 0) { 
k++; 

return 1; 
} else { 
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fksinf(&fl, &f2, &f3, k) ; 

ret = fl*ak_l + f2*ak_2 + f3*ak_3; 

} 

ak_3 = ak_2; 
ak_2 = ak_l; 
ak_l = ret; 

k++; 

return ret; 



// Scalar wave, zero frequency 

define D(a) double (a) 

typedef complex<double> cmplx ; 

cmplx Ylm(int 1, int m, double theta, double phi); 

void test_zero (double radius); 

void getZeroPsiH (cmplx* PsiH, cmplx* dPsiH) ; 
void getZeroPsil (cmplx* Psil, cmplx* dPsil) ; 

const double PI(M_PI); 

using namespace ScalarWave ; 



double Bang2LDL2(int L) { 
assert (L >= 0) ; 
double ret = 1 ; 
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for (int i = L; i > 0; i~) { 
ret *= (2.*i)*(2.*i-l)/(i*i) ; 

} 

return ret; 



double AnHorizon (double* Aninl, int n, int 1) { 
return *Aninl = -(*Aninl/(n*n))*(l+n)*(l-n+l) ; 

} 



double Aninfinity (double* Anml, int n, int 1) { 
*Anml = (*Aninl/(n*(2*l+l+n)))*(l.*(l+n)*(l+n)) ; 
return *Aninl; 

} 



// Wronskian/M for zero frequency scalar modes 

double WzeroDM(int 1) { 

return (iy„2==0?-l:+l)*Bang2LDL2(l)*2* (2*1+1) ; 

} 



void getZeroPsiH(cmplx* PsiH, cmplx* dPsiH) { 

*PsiH = 0; 
*dPsiH = 0; 

double b = 1 ; 

for (int n = 0; n <= 1; n++) { 
cerr « SHOW(b) « endl; 
double Rn = pow(R/(2*M) ,n) ; 
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*PsiH += b*Rn; 
*dPsiH += n*b*Rn/R; 

cout « SHOW(*PsiH) « SHOW(*dPsiH) « endl; 
AnHorizon(&b,n+l, 1); 



assert (b == 0) ; 

cout « SHOW(*PsiH) « SHOW(*dPsiH) « endl; 



void getZeroPsiI(cmplx* Psil, cmplx* dPsil) { 

*PsiI = 0; 
*dPsiI = 0; 

double a = 1 ; 

for (int n = 0; n <= 100; n++) { 
cerr « SHOW(a) « endl; 
double Rn = pow(2*M/R,n) ; 
*PsiI += a*Rn; 
*dPsiI += -n*a*Rn/R; 

cout « SHOW(*PsiI) « SHOW(*dPsiI) « endl; 
if (n>l) { 

double corr = abs(a*Rn/(*PsiI)) + abs(-n*a*Rn*R/(2*M)/(*dPsiI)) ; 
if (corr < l.e-17) { 

cerr « SHOW(n) « SHOW(corr) « endl; 

break; 

} 

} 

Aninf inity(&a,n+l,l) ; 



*PsiI *= pow(R/(2*M) ,-(1+1)) ; 

*dPsiI = pow(R/(2*M),-(l+l))*(*dPsiI) -(l+l)*(*PsiI)/R; 



128 



} 



void test_zero() { 



cout << "Entering test_zero" « endl; 
time_t tenter = time(O); 



// Set ScalarWave variables: 



M = 1; 
m = 0; 



// mass of the hole 

// spherical harmonic index 



R = 20*M; // radius of the orbit 

w = M/ (R*R*R) ; // angular frequency of orbit 



for (int 1=0; 1 < 41; 1++) { 

cmplx PsiH; 
cmplx dPsiH; 

getZeroPsiH(&PsiH, fedPsiH) ; 

cmplx Psil = 0; 
cmplx dPsil = 0; 

getZeroPsil (fePsil , fedPsil) ; 

cout « "Here is a Wronskian test" << endl; 

cmplx W = R*R*(PsiH*dPsil - PsiI*dPsiH)*(l-2*M/R) ; 

cout « setprecision(12)« SHOW(l) « SHOW(W) 

« SHOW(M*WzeroDM(l)) « endl; 
cout « SHOW(l) « SHOW(abs((W - M*WzeroDM(l) ) /W) ) 

« SHOWCW - M*WzeroDM(l)) « endl « endl; 

cout « SHOW(Psil) « SHOW(dPsil) « endl; 
cout « SHOW(PsiH) « SHOW(dPsiH) « endl; 



omega = m*w; 



// frequency of the radiation 



} 
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tiine_t texit = time(O); 

cout« "test_zero execution time is " « diff time (texit, tenter) 
<< " seconds" « endl; 

cout « "Leaving test_zero" « endl; 
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